探討石墨板與奈米碳管之機械性質

國 立 交 通 大 學

機械工程研究所

碩 士 論 文

探討石墨板與奈米碳管之機械性質

Investigating Mechanical properties of Graphene sheet and

Carbon Nanotubes

研 究 生 ：涂潔鳳

指導教授 ：蔡佳霖 博士

探討石墨板與奈米碳管之機械性質

Investigating Mechanical properties of Graphene sheet and

Carbon Nanotubes

研 究 生：涂潔鳳 Student：Jie-Feng Tu 指導教授：蔡佳霖 Advisor：Jia-Lin Tsai

國 立 交 通 大 學

機 械 工 程 研 究 所

碩 士 論 文

A Thesis

Submitted to Department of Mechanical Engineering College of Engineering National Chiao Tung University

in partial Fulfillment of the Requirements for the Degree of Master

in

Mechanical Engineering

June 2009

Hsinchu, Taiwan, Republic of China

探討石墨板與奈米碳管之機械性質

研究生：涂潔鳳 指導教授：蔡佳霖 博士

國立交通大學機械工程學系碩士班

摘 要

本 研 究 利 用 分 子 動 力 學 理 論 分 析( Molecular Mechanics Analytical Solution )與分子動力學( Molecular Dynamics )探討石墨板與奈米碳管的機 械性質，碳-碳原子間的力場則引用 AMBER 勢能[1]描述之。在石墨板分析 中，比較解析解與分子動力學兩種方法分析的結果。另外在分子動力學的 模擬中，將石墨板施以平面方向的應力或應變來求得楊氏模數( Young’s modulus )以及浦松比(Poisson ratio)。分析結果顯示，利用解析解與分子動力 學所求得之平面方向材料模數，均符合 G=E/(1+2υ)等向材料之特性。出平 面方向機械性質主要利用Lennard - Jones 勢能[2]來探討剪力模數，對石墨 板施以出平面方向剪應力，產生層與層間相對位移，而得到出平面方向剪 力模數。比較平面方向與出平面方向之機械性質，發現出平面方向之機械 性質均比平面方向機械性質要低許多。最後討論不同文獻所引用凡得瓦勢 能的參數對於石墨板理論結果的影響。 奈米碳管的研究中，則利用分子動力學模擬並分析單壁碳管和雙壁碳 管的楊氏模數及浦松比，研究鋸齒形(zigzag)與扶手椅形(armchair)兩種不同 原子排列結構，以及不同管徑下對於機械性質之影響。發現在管徑較小時， 鋸齒型(zigzag)碳管的楊氏模數會比扶手椅形(armchair)碳管要來的低；而在 管徑較大時，則不會受到原子結構排列的影響。接著利用雙壁碳管內管與

外管相互滑動時，所產生的層間原子作用力(atomstic interlayer force)，模擬 方法分為兩種，(1) 將雙壁奈米碳管的內管拉出，並計算內管所需之合力。 (2) 利用兩種不同平衡條件的雙壁碳管，使內外管相對旋轉後，計算碳管的 內管在旋轉時所需之扭轉力，再觀察扭轉力與能量的變化。並藉由層間作 用力的分析，解釋在奈米複合材料的應用中，雙壁奈米碳管層間應力傳遞 的能力。

Investigating Mechanical properties of Graphene sheet and

Carbon Nanotubes

Student：Jie-Feng Tu Adviser：Dr. Jia-Lin Tsai

Dependent of Mechanical Engineering

National Chiao Tung University

Abstract

This study aims to investigate the mechanical properties of graphene and single/multi-walled carbon nanotubes(S/MWCNT) by analytical molecular structural mechanics and molecular dynamics (MD). The interatomic covalent bond of carbon-carbon is represented by AMBER force field which was employed in the analytical molecular mechanics and MD simulations as well, and the non-bonded interaction was described by Lennard-Jones potential. In particular, in-plane elastic properties of MD graphene model obtained from applying stress or strain were compared with analytical solutions. Results indicated that in-plane properties, predicted by analytical solution demonstrate a good agreement with MD solution. It was found that in-plane properties of graphene basically

satisfied the relation, G₁₂ = E₁/(1+2ν), and thus it can be regarded as a transverse isotropic

material. Applying stress on the graphene model in the out-of-plane direction and calculating the relative displacement between graphene layers interacted by van der Waals potential, out-of-plane shear modulus could be derived.

Because of the weak Lennard-Jones potential, the in-plane mechanical properties are higher than out-of-plane shear modulus. Since van der Waals potential may affect the mechanical properties of the graphene model, the effect of van der Waals on graphene mechanical properties was also investigated in the study.

The axial Young’s modulus and Poisson ratio of S/DWCNTs is predicted by the MD method. In this article, elastic properties of unchiral (zigzag and armchair) S/DWCNTs of different diameters are studied and discussed. The results indicated that Young’s modulus of

CNTs varies with the tube diameter and is affected by their helicity. With increasing the tube diameter, the Young’s modulus of both armchair and zigzag CNTs increases. The Young’s modulus of zigzag CNTs is lower than that of armchair at smaller diameters, however, at larger diameters there is no significant difference observed between armchair and Zigzag CNTs.

The interlayer atomistic force of between the neighboring graphene layers in DWCNTs was examined by performing inner tube pullout and rotation tests with respective to the outer tube at stationary position. The pullout force corresponding to the pullout distance was measured in the pullout test and the torque of inner tube and interlayer VDW energy variation associated with the rotation angle was calculated in the rotational tests. This information is very essential for understanding the load transfer efficiency between the graphene layers in DWCNTs when they was considered an reinforcement in the nanocomposites.

誌 謝

經過研究所的磨練，讓我獲益良多，也深深地體驗到認真、以及負責 任的態度是非常重要的。在此我要感謝指導教授 蔡佳霖博士，在他細心 以及耐心的指導下，使我不僅順利完成論文的研究，更教給我正確且積極 的態度。在此，同時也感謝清華大學動機械葉孟考教授、交通大學機械系 蕭國模教授撥空擔任學生口試委員，並且在口試的過程中給予保貴之意 見，讓此篇論文更佳完整。 此外要感謝實驗室的學長曾世華、齊揚楷、蕭宏給予研究上的建議和 指導，以及研究室所有的學弟以及同學，有你們讓我的研究生活增添許多 樂趣。 感謝我的家人以及好友給予我的支持，尤其是兩位妹妹平時的關心和 鼓勵，讓我更放心的去追求自己的目標。

目 錄

中文摘要 ……… i 英文摘要 ……… iii 誌謝 ……… v 目錄 _{……… Vi} 表目錄 ……… _ix 圖目錄 ……… _x 第一章 緒論………. ₁ 1.1 研究背景與文獻回顧………..………... ₁ 1.2 研究目標….………. 11 第二章 石墨板………. 12 2.1 石墨板結構……….……… 12 2.2 原子間之勢能函數……….………. 12 2.3 分子力學理論分析………. 13 _{2.3.1 平面載重分析……….…….……..………. 13 2.3.1.1 平面單軸方向拉伸載重：楊氏模數、浦松比 15 2.3.1.2 平面剪力載重:平面方向剪力模數………….. 17 2.3.2 出平面方向載重分析:出平面方向剪力模數……….. 19} 2.4 分子動力學模擬………...………... 25

_{2.4.1 分子動力學(MD)模擬……… 25} 2.4.2 石墨板模型………. _{29 2.4.3 石墨板平面方向機械性質………. 32 2.4.3.1 平面單軸方向拉伸載重：楊氏模數、浦松比 32 2.4.3.2 平面剪力載重:平面方向剪力模數………….. 34 2.4.4 出平面方向載重分析:出平面方向剪力模數………... 35} 2.5 結果與討論……….. 36 _{2.5.1 石墨板材料模數與參考文獻之比較與討論……….… 36 2.5.2 不同凡得瓦力勢能對於模擬之結果及影響…………. 38} 第三章 奈米碳管……….………... 40 3.1 奈米碳管結構及勢能……….. 40 3.2 奈米碳管軸方向機械性質………... 42 3.2.1 單壁奈米碳管(SWCNT)MD 模型………. 42 3.2.2 單壁奈米碳管軸方向機械性質…..………. 43 3.2.3 雙壁奈米碳管 MD 模型…..……….………. 45 _{3.2.4 雙壁奈米碳管軸方向機械性質.….………. 46} 3.3 雙壁奈米碳管層間原子作用力……….……… 48 3.3.1 有限長度雙壁奈米碳管 MD 模型………. 48 3.3.2 雙壁碳管層間原子作用力…..………. 49

3.4 結果與討論……….……….……… 52 _{3.4.1 奈米碳管之楊氏模數比較………. 52} 3.4.2 雙壁碳管層間原子作用力…..………. 54 第四章 結論與未來展望……… 57 4.1 結論……….……… 57 4.2 未來展望……….……….……… 61 參考文獻 _{…….……… 62}

表 目 錄

表 _{2.1 石 墨 板 模 型 (2) 、 (3) 受 到 單 方 向 應 變 求 得 之 勁 度 矩 陣 係} 數……….……. 66 表 2.2 平面方向材料參數(E1、υ12)………...…... 66 表 _{2.3 平面方向剪力係數材料參數(G12)……….……… 67} 表 _{2.4 出平面方向剪力係數材料參數(G13)………. 67} 表 2.5 不同凡得瓦力勢能之參數比較表………. 68 表 _{2.6 利用不同凡得瓦力勢能之參數在模型(2)所得之平面方向材料參} 數(E1、υ12、G₁₂)………...…. 69 表 _{2.7 利用不同凡得瓦力勢能之參數在模型(3)所得之平面方向材料參} 數(E1、υ12、G₁₂)….………. 70 表 2.8 利用不同凡得瓦力勢能之參數在模型(3)所得之出平面方向材料 參數(G13)……….. 71 表 3.1 各種不同半徑之碳管模型尺寸………...……….. 72 表 3.2 利用各種條件下單壁奈米碳管模型 A、B 所得之軸方向材料參 數……….. 73 表 3.3 (a) (5,0)(14,0)雙壁奈米碳管模型與單壁奈米碳管模型軸方向材 料參數之比較……….. 74 表 3.3 (b) (3,3)(8,8)雙壁奈米碳管模型與單壁奈米碳管模型軸方向材料 參數之比較……….. 74

圖 目 錄

圖 _{2.1 石墨板結構………. 75} 圖 _{2.2 石墨板 A-B 疊層結構………..………...…... 75} 圖 2.3 原子間作用力的形式……….……… 76 圖 2.4 (a) 單方向拉伸載重 (b) 剪力載重……….………. 76 圖 2.5 簡化石墨模型……….……….... 77 圖 2.6 簡化石墨模型受拉伸均佈載重作用等效為集中載重作用…...….. 77 圖 2.7 OB 鍵結長度與角度變化幾何關係變化圖…..………. 78 圖 2.8 石墨模型受到均佈剪力載重等效為集中載重………..….. 78 圖 2.9 鍵結長度與角度變化幾何關係變化圖………….……...……...….. 79 圖 2.10 單顆原子對於鄰近層之截斷半徑示意圖……...……….. 79 圖 2.11 _{碳原子之距離對應凡得瓦勢能與凡得瓦力…...……….. 80} 圖 _{2.12 單顆碳原子相對於相鄰層分為}

α

_與

β

_{兩種平衡位置…...……….. 80} 圖 _{2.13 原子位置改變關係………...……….. 80} 圖 _{2.14 單層石墨板不考慮凡得瓦力(模型 1)……….……….. 81} 圖 2.15 單層石墨板且考慮凡得瓦力(模型 2)……….……….. 81 圖 2.16 八層石墨板考慮凡得瓦力(模型 3)…….……….. 82 圖 2.17 石墨板(模型 1)不受力之(a)能量變化圖、(b)溫度變化圖、(c)1 方 向壓力變化圖、(d)2 方向壓力變化圖、(e)3 方向壓力變化 圖……….……….………….. 83 圖 2.18 石墨板(模型 2)不受力之(a)能量變化圖、(b)溫度變化圖、(c)1 方 向壓力變化圖、(d)2 方向壓力變化圖、(e)3 方向壓力變化 圖……….……….………….. 84 圖 2.19 分子動力學模擬之石墨板(模型 1)………...……….. 85

圖 2.20 MD 石墨結構(模型 1-3)施以單軸拉伸應力……….. 85 圖 _{2.21 石墨板(模型 1)受單軸應力之(a)能量變化圖、(b)溫度變化圖、(c)1} 方向壓力變化圖、(d)2 方向壓力變化圖、(e)3 方向壓力變化 圖……….……….………….. 86 圖 2.22 MD 石墨結構(模型 2-3)施以單軸拉伸應變.……….……….. 87 圖 2.23 MD 單層石墨(模型 1)施結構施以剪應力…….…….……….. 87 圖 2.24 MD 單層石墨(模型 2-3)施結構施以剪應力…..……….. 87 圖 2.25 MD 多層石墨結構(模型 3)施以出平面剪應力……..……….. 88 圖 3.1 以向量表示石墨捲曲方式說明圖[30]…...………….……….. 89 圖 3.2 勢能函數之曲線擬合(curve fitting)曲線圖………….……….. 89 圖 _{3.3 奈米碳管之雙面角結構圖……….……….……….. 90} 圖 _{3.4 單壁奈米碳管模擬室……….……….……….. 91} 圖 _{3.5 單壁奈米碳管(5,0)不受力狀態下之(a)能量變化圖、(b)溫度變化} 圖、(c)1 方向壓力變化圖、(d)2 方向壓力變化圖、(e)3 方向壓力 變化圖……….……….……….. 92 圖 3.6 單壁奈米碳管(5,0)在拉伸方向(Z 方向)對模擬室施以應力 σBox 之(a)能量變化圖、(b)溫度變化圖、(c)1 方向壓力變化圖、(d)2 方向壓力變化圖、(e)3 方向壓力變化圖………... 93 圖 3.7 SWCNT 在原子排列結構不同時單壁碳管對於管徑大小改變之 楊氏模數……….………….……….. 94 圖 3.8 雙壁奈米碳管(5,0)(14,0)分子動力模擬下之結構.……….. 95 圖 3.9 雙壁奈米碳管(5,0)(14,0)在不受力狀態下之(a)能量變化圖、(b) 溫度變化圖、(c)1 方向壓力變化圖、(d)2 方向壓力變化圖、(e)3 方向壓力變化圖……….…….…….……….. 96 圖 3.10 雙壁奈米碳管(5,0)(14,0)在拉伸方向(Z 方向)對模擬室施以應 力σBox之(a)能量變化圖、(b)溫度變化圖、(c)1 方向壓力變化圖、 (d)2 方向壓力變化圖、(e)3 方向壓力變化圖……….……….. 97

圖 3.11 1 軸向非連續性之雙壁碳管(a) (5,0)(14,0) (b) (9,9)(14,14)...…….. 98 圖 _{3.12 雙壁碳管結構未受束縛的平衡狀態.……….………... 99} 圖 _{3.13 非週期性雙壁碳管之內管沿軸方向平移拉出結構圖...…….. 100} 圖 _{3.14 非週期性雙壁碳管之拉出碳管之力對應拉出平移長度……..….. 100} 圖 _{3.15 雙壁碳管內管逐漸拉出觀察圖……..….……...……….. 101} 圖 3.16 雙壁碳管拉出內管時之內外管結構之相對位置…………..……. 102 圖 3.17 雙壁碳管做內管剛體旋轉之模型……..……….. 103 圖 3.18 雙壁碳管(5,0) (14,0)在平衡條件(A)下之(a) 扭轉力變化圖 (b) 層間凡得瓦勢能變化圖……..………..….... 104 圖 3.18 雙壁碳管(5,0) (14,0)在平衡條件(A)下之 (c) 內管中心軸 (d)外 管中心軸………..……….. 105 圖 _{3.19 雙壁碳管 (5,0) (14,0)在平衡條件(B)下之(a) 扭轉力變化圖 (b)} 層間凡得瓦勢能變化圖……..………..….... 106 圖 3.19 雙壁碳管(5,0) (14,0)在平衡條件(B)下之 (c) 內管中心軸 (d)外 管中心軸………..……….. 107 圖 3.20 雙壁碳管(9,9) (14,14)在平衡條件(A)下之(a) 扭轉力變化圖 (b) 層間凡得瓦勢能變化圖……..………..….... 108 圖 3.20 雙壁碳管(9,9) (14,14)在平衡條件(A)下之 (c) 內管中心軸 (d) 外管中心軸………..……….. 109

第一章 緒論

1.1 研究背景與文獻回顧

近幾十年來，複合材料不論是理論研究或者實際量產化商品，不斷推 陳出新，廣泛受到應用。複合材料主要由基材(Matrix)及補強材(reinforcement) 建構而成，在傳統商用複合材料中使用的補強材尺寸，大都為半徑在數微 米(μm)至數百微米之間，其中最被廣泛使用的是碳纖維(carbon fiber)。而目 前炙手可熱的奈米材料，如奈米碳管(Carbon Nanotube, CNT)等，補強效果 比以往的碳纖更好，使它成為具有潛力的先進複合材料，另外奈米材料的 光、電、磁及其他物理化學的性質大多會異於巨觀尺度的材料性質，經過 十幾年的發展，奈米材料應用在複合材料上的發展潛力已逐漸浮現。 以機械性質而言，奈米碳管是一個非常強韌的物質，它的機械強度非 常好而且具有可反覆彎曲後並不斷裂的特性，以單壁奈米碳管為例，它的 強度約為鋼的10-100 倍，但是重量卻只有鋼的1/6，是一種輕且機械強度非 常好的材料。 奈米碳管以及石墨板中基本為sp2_{混成形成的碳-碳共價鍵，其楊氏模數} 與sp2 共價鍵密切相關，所以可以得知原子與原子之間的關係對機械性質有 一定之影響。 在過去的研究文獻中，碳管的發現首先由Iijima等人[3]經由碳電弧放電 法(Carbon arc discharge)，發現了一些針狀物，利用高解析穿透式電子顯微

鏡(High Resolution Transmission Electron Microscope)觀察後，發現針狀物為 多層同軸中空管，主要由碳原子組成，此結構即被稱為多壁奈米碳管，亦 觀察出多壁碳管在層與層之間之間距為3.4Å。理論方面。之後亦有學者[4-7] 針對石墨以及奈米碳管的機械性質的理論解析解推導有相當的貢獻，Xiao 等人[4]利用Morse勢能配合理論解方法來預測單壁碳管方面之軸方向楊氏 模 數(Young’s modulus) 、 浦 松 比 (Poisson ratio) 以 及 剪 力 模 數 (Shear modulus)。將碳管原子模型取出代表性的單元，施以固定應變，並求得鍵結 伸長量以及鍵結角度變化之關係式；利用力平衡推導出鍵結伸縮力與鍵結 角度彎矩之關係； 而將勢能對變形量微分，可得到力與變形量的關係式；結合各種關係式， 則可得到外力與原子模型變形之關係，即得到碳管的材料係數。由理論解 發現碳管在曲率越小時，則楊氏模數會越大；並且在碳管直徑小於12Å時， 扶手椅形(armchair)單壁碳管的楊氏模數要比鋸齒形(zigzag)單壁碳管明顯 高出許多；而管徑大於12Å時，扶手椅形(armchair)單壁碳管與鋸齒形(zigzag) 單壁碳管的楊氏模數幾乎相同。Kelly等人[5]以兩層石墨板為石墨模型，利 用解析解求得石墨板的出平面材料模數，並將原子與相鄰石墨板上的原子 分為兩種相對位置，引用最小勢能原理，引用Lennard-Jones勢能計算出能量 密度(energy density)對應變微分，可得到石墨板的勁度矩陣參數C33、C44， 與實驗結果做比較，可獲得相當接近的結果。Cho等人[6]沿用前面敘述之解

析解方法[4-5]，並利用Morse勢能以及二維碳原子模型，假設在微小變形量 下，推導出石墨板平面方向的楊氏模數(in-plane Young’s modulus)、浦松比 (in-plane Poisson ratio)和剪力模數(in-plane Shear modulus)，其理論解結果符 合等向性材料G₁₂ = E₁/2

(

1+υ₁₂

)

之定則，並與實驗值比較可得到相當好的比對 結果。而在出平面方向材料模數的推導則利用最小勢能原理，將能量密度 對應變微分，可得到勁度矩陣中的C₃₃、C₄₄與出平面方向浦松比υ₁₃，再進而 推算出平面方向楊氏模數E₃、剪力模數G₁₃。結果顯示石墨板在平面方向的 強度會比出平面方向增加許多，而υ₁₃近乎於零，表示石墨板受到平面方向 的載重時，出平面方向變形很小，同樣地，石墨板受到出平面方向載重時， 亦對平面方向應變影響不大。Scarpa [7]等人則利用桁架模型以及蜂巢模型 建立石墨板的代表性元素，並利用AMBER與Morse兩種不同的勢能參數分 析石墨板的平面方向機械性質，結果顯示引用AMBER模擬石墨板平面方向 的材料模數會比Morse力場的結果低。而Reddy等人[8] 引用Tersoff-Brenner 勢能，利用最小能量平衡的方法，將有限長度之石墨板在應力為零的穩定 狀態下給予位移，得到石墨板在不同位移下的能量，再引用Cauchy-Born rule 得到石墨板之材料模數。實驗方面，Blakslee等人[9]利用超音波波動原理， 用探頭給予石墨試片擾動測得石墨板在各方向之波速，以得到石墨板之材 料模數。Al-Jishi等人[10]利用聲子(Phonon)所釋放之能量測得石墨板的各種 模態，透過材料特性之公式解得到石墨板之材料模數。Dolling及Jenkins等

人[11-12]則透過不同的實驗模型，測量石墨板之震動頻率來以求得其材料 特性。文獻中使用解析解和分子動力學模擬，所求得的石墨材料係數均接 近實驗的結果。 目前已有許多學者[13-19]利用分子動力學模擬石墨板或奈米碳管的材 料模數，並且接近實驗結果。部分文獻中，假設單壁碳管壁厚約3.4Å [13, 16, 17, 19, ]，並且探討管徑大小以及原子排列結構對於分析結果之影響。Bao 等人[13]利用分子動力學方法，研究石墨板和單壁碳管在不同原子排列結構 下的楊氏模數，藉由給定一外力，利用REBO(reactive empirical bond-order) 勢能以及Lennard-Jones勢能計算碳管軸向伸長變化量，以應力應變關係得到 軸方向之楊氏模數(Young’s modulus)。結果顯示，石墨板在不同層數(1~5 層)、不同面積的條件下，會得到不同的楊氏模數，值約在1012~1033GPa之 間，平均值為1026GPa與實驗值1025GPa相近，表示理論值已成功地接近實 驗值；而碳管在不同長度下，楊氏模數的理論值介於929.8±11.5GPa，比石 墨板稍微小一些。另一方面，單壁碳管在不同的原子排列結構下，理論結 果顯示扶手椅形(armchair)、鋸齒形(zigzag)、手性形(chiral)碳管的楊氏模數 均有相似的理論值，表示原子排列結構和管徑大小不會特別影響碳管機械 性質。Kohji等人[14]引用Tersoff-Brenner 勢能計算，同樣假設碳管在微小 變形下，利用兩種方式計算碳管楊氏模數。第一種方式為施以應變法，將 碳管軸方向模擬為周期性邊界，表示碳管軸被想像為無限長之模型，將應

變分為N個步驟施加，利用共軛梯度法計算(the conjugated gradient method) 收斂之能量最小值，並利用最小能量原理，求得碳管直徑改變量；另一種 方式是在有限長度的碳管上施軸向力，利用能量平衡，得到碳管伸長量， 兩種方法均可得到合理的碳管楊氏模數(Young’s modulus)和浦松比(Poisson ratio)。且假設兩種不同碳管截面的面積公式，用來模擬單一碳管之材料性 質以及複合材料中添加之奈米碳管的楊氏模數。結果顯示，在單壁碳管管 徑約大於5Å時，其楊氏模數則會越大，直到趨近於相同的值；而原子排列 結構對於楊氏模數的影響方面，扶手椅形(armchair)的碳管的楊氏模數在管 徑較大時比鋸齒形(zigzag)高出一些。Zhang 等人[15]模擬單壁碳管受到壓 力時的狀態，且利用挫屈理論公式推算出楊氏模數並得到管壁厚度，結果 顯示單壁碳管壁厚比一般的3.4Å小，利用此方法所推算出之楊氏模數會較 高。Agrawal 等人[16]用分子動力學模擬透過五種方法來求得單壁碳管的楊 氏模數和浦松比，並研究凡得瓦力對於碳管機械性質之影響，結果顯示利 用側向振動頻率所求出之楊氏係數值會與其他方法之結果差異較大，而利 用應力法、應變法、應變能法、軸向振動頻率測法所得到的(14,14)單壁碳 管楊氏模數大約在0.73~0.76TPa，而凡得瓦力對於單壁碳管機械性質幾乎沒 有影響。並模擬直徑介於9.5~19Å的扶手椅形(armchair)及鋸齒形(zigzag)單 壁碳管，結果顯示相同管徑下，扶手椅形(armchair)會比鋸齒形(zigzag)單壁 碳管的楊氏模數要高出一些，且管徑越大時，其楊氏模數則會越小，最後

會趨於一定值。Chen等人[17]利用分子動力學探討單壁碳管在不同原子排列 結構下，考慮凡得瓦力(van der Waal atomistic interaction)與忽略凡得瓦力間 的差異，結果證明在考慮凡得瓦力時可使楊氏模數(Young’s modulus)增加 9%，浦松比(Poisson ratio)降低27%，剪力模數(Shear modulus)上升12%，且 在多壁碳管中，碳管內層剪力強度增加13%。證明了凡得瓦力在分子動力學 上對於模擬碳管機械性質的影響是一項值得注目的研究；而原子排列結構 的影響方面，結果顯示在相同管徑下，扶手椅形(armchair)碳管會比鋸齒形 (zigzag)單壁碳管的楊氏模數要高出一些，而管徑越大時，其楊氏模數則會 越大，最後會趨於一定值。而在利用其他方法分析碳管機械性質的文獻中， Natsuki 等人[18]則結合分子力學及固體力學成為一理論解析方法，來求解 單壁碳管在不同原子結構排列以及管徑下之楊氏模數、浦松比和剪力模 數，且此文獻的解析解結果顯示，在相同管徑下，原子結構排列不影響碳 管機械性質，而碳管管徑越大時，其楊氏模數會越低，直到趨於一穩定值。 Li 等人[19]則將原子間之鍵結以梁元素取代，建立單層石墨或單壁碳管結 構，求得勁度矩陣，並施以單軸力來得到楊氏模數、浦松比、剪力模數， 結果顯示單壁碳管管徑越大時，碳管的楊氏模數越大且接近石墨板；在管 徑小於約7Å時，扶手椅形(armchair)碳管的楊氏模數會比鋸齒形(zigzag)要高 一些，而管徑大於約7Å時則相反。Tserpes 等人[20]利用有限元素法建立單 壁碳管模型，給定單軸拉力或是扭轉力，並比較碳原子在不同結構排列下

之模擬結果，其結果顯示楊氏模數會受到原子排列結構和管徑大小之影 響，並在管徑越大則楊氏模數則越大；而在原子排列結構的影響方面，管 徑小於10Å的單壁碳管，扶手椅形(armchair)碳管的楊氏模數會比鋸齒形 (zigzag)要高一些，管徑較大時則相反。Lu [21]利用理論模型公式以能量對 應變作二次微分求得單壁碳管以及石墨板之材料參數，單壁碳管的分析結 果顯示，楊氏模數不受原子排列結構及管徑大小的影響。整理上述文獻， 探討原子排列結構不同時，對於楊氏模數的影響，發現部分文獻[4, 19, 20] 單壁碳管在管徑較小時，扶手椅形(armchair)單壁碳管的楊氏模數比鋸齒形 (zigzag)單壁碳管要高，但是管徑較大時則發現扶手椅形單壁碳管的楊氏模 數比鋸齒形單壁碳管低；另外，亦有文獻[14, 16]發現單壁碳管不論管徑大 小，鋸齒形(zigzag)碳管的楊氏模數一律比扶手椅形(armchair)碳管要高；其 他部分文獻[13, 15, 18, 21]表示在相同的管徑下，碳管的楊氏模數不會受到 原子結構排列的影響。而在討論管徑大小對於楊氏模數的影響方面，有部 分文獻[4, 14, 17, 19, 20]研究顯示，在管徑越大時，其楊氏模數則越大；另 外亦有文獻[16, 18]指出，在管徑越大時，單壁碳管的楊氏模數則越低；其 他部分文獻[13, 21]指出，管徑大小並不影響碳管之楊氏模數；上述研究單 壁碳管的文獻全部指出，楊氏模數在管徑曲率半徑越大時，其楊氏模數則 會越接近石墨板；。 多壁碳管的研究方面，亦有多位學者探討其機械性質[17, 21, 23]，以及

研究雙壁碳管的外管單獨受到拉伸時，凡得瓦力傳送拉伸力量至內管的能 力[22, 23]。Chen等人[17]針對雙壁碳管在不同原子排列結構和有無凡得瓦 力的情形下，研究楊氏模數所受到的影響，結果顯示雙壁碳管在考慮凡得 瓦力時，其楊氏模數則會提高；且不論有無考慮凡得瓦力時，扶手椅形 (armchair)雙壁碳管的楊氏模數會比鋸齒形(zigzag)雙壁單壁碳管要高一 些，而手性形(charial)碳管的楊氏模數則介於扶手椅形(armchair)和鋸齒形 (zigzag)兩者之間；此篇文獻結果顯示，單壁碳管和雙壁碳管有相同的趨勢。 Lu [21]利用理論模型公式模擬雙壁碳管，而此文獻指出單壁碳管和雙壁碳 管的楊氏模數，均不受管徑大小和原子排列結構的影響。Shen 等人[22]利 用分子動力學研究單軸力單獨施於外管碳管時，是否可以透過凡得瓦力的 傳送而使內管同樣有受到拉伸或壓縮之效果，結果顯示當單軸應變僅施於 外管時，內管所受到的力卻非常微小，表示凡得瓦力對於內外管之間拉伸 的傳送能力非常的有限。Li 等人[23]引用前面文獻[19]建立單壁碳管之模型 之方法建立雙壁碳管，而原子間的凡得瓦力則使用桁架模擬之。並研究多 壁碳管及多層石墨板的楊氏模數、剪力模數，以及雙壁碳管內管與外管層 間凡得瓦力傳遞力量之能力，結果顯示，管徑相同的雙壁碳管和單壁碳管， 兩者楊氏模數相當接近，可見雙壁碳管的機械性質，會受到所組成的單壁 碳管的影響；當單軸力僅只有施於外管時，透過凡得瓦力傳送到內管的力 非常微小，表示凡得瓦力對於內外管之間的作用力非常微小。整理上述關

於雙壁碳管文獻的結論，可以看出雙壁碳管與單壁碳管的楊氏模數有很大 的關聯；而凡得瓦力對於雙壁碳管層間力量的傳遞，則扮演較薄弱的角色， 由此可觀察出凡得瓦力對於雙壁碳管楊氏模數的影響不大。 雙壁碳管層間作用力的研究中，有文獻[17, 24-26]討論在拉出或轉動內 管時，內外管相對位置與層間作用力之關係，Chen等人[17]探討不同原子排 列結構下，內管拉出距離與其拉出內管所需合力之關係，結果顯示在相同 管徑下，不同的原子排列結構對於拉出內管之力並沒有影響，但內管拉出 力的曲線會因原子排列結構不同而呈現不同的震盪情形，其振幅大小亦受 到原子排列結構的影響，震盪曲線的平均值顯示拉出碳管時，一開始會逐 漸增大，拉出一固定距離之後，則會呈現一定值，直到內管完全拉出時則 降為零。Tserpes [24]研究利用有限元素法建立雙壁碳管模型，控制內層與 外層間距為3.4~5.0Å之間，並計算拉出內管碳管時所需之力。結果顯示當內 外層間距越大時，其所需之力則越小，且同樣在抽出碳管時，拉出之力一 開始會呈現線性現象，最後則為一定值。Xia等人[25]利用分子動力學研究 雙壁碳管的機械性質，且內管碳管分為末端為封閉及非封閉兩種，在拉出 內管碳管時，層與層間相對滑動所產生原子層間作用力即視為磨擦力，亦 看作內管碳管移動時，內外管層間所產生之剪力，並且利用分析結果解釋 實驗之現象。碳管管壁上分佈的均佈壓力，會使內管滑動時所產生更大的 摩擦力，讓抽出碳管之力增加，而內部碳管逐漸脫離碳管時，抽出碳管之

力則減少。Zhang 等人[26]利用分子動力學研究雙壁碳管在轉動內管時的最 小勢能變化量，同時探討高溫時，轉動碳管對於能量變化的影響，結果顯 示在轉動內管碳管時，最小勢能呈現穩定的週期性變化，並可利用公式計 算其能量週期長度。轉動內管時的能量變化量，不論在最小勢能狀態或是 升溫狀態時，相較於拉出碳管之力均非常的微小。由上述文獻可知，雙壁 碳管結構的原子層間作用力，透過原子間凡得瓦力的影響而產生，且對於 機械性質的剛性影響非常微小。而在拉出內管碳管的過程中，在一開始拉 碳管之力會逐漸增加，直到拉出一固定距離，拉出內管碳管之力則會呈現 一定值，最後當內管完全拉出時，拉力會降為零，表示內外管幾乎沒有凡 得瓦力的作用。

1.2 研究目標

本研究目的為模擬出奈米石墨板及奈米碳管之機械性質。第二章引用 Amber勢能[1]之參數，主要利用Xiao[4]理論堆導方法，將能量對變形量微 分而得到鍵結上的作用力，利用力平衡與幾何變形關係推求石墨板在平面 方向之材料參數E1、υ12、G12。參考Kelly[5]的推導方法，以Lennard-Jones 勢能計算能量密度並對應變微分，而得到C44(=G13)。同樣地，利用分子動力 學理論，引用Berendsen[27]之數值方法，以及Amber勢能，在石墨板分別施 以單方向應力，而得到石墨板所改變之1、2方向應變，利用應力應變關係 求得石墨板之材料參數E₁、υ₁₂、G₁₂、G₁₃。另一種方法則利用勁度矩陣法， 將石墨板施以一應變，求得勁度矩陣中之參數，再求得石墨板之楊氏模數。 最後探討解析解與分子動力學兩種方法的差異。 第三章則利用分子動力學模擬分析單壁碳管的楊氏模數及浦松比，並 研究鋸齒形(zigzag)與扶手椅形(armchair)兩種不同原子排列結構，以及不同 管徑對於碳管機械性質的影響。接著建立雙壁碳管模型，探討其管壁之間 的層間原子作用力，模擬方法是將雙壁奈米碳管做內管拉出或內外管相對 旋轉，最後得到內外碳管相對位移以及內管所受之力的關係。 最後第四章將第二、三章的分析結果做比較，並且探討其機械性質。

第二章 石墨板

2.1 石墨板結構

石墨板(Graphene)主要是以碳原子組成，其碳-碳原子間鍵結會有 3 個 sp2軌域與一個 p 軌域，可在碳-碳間各取一 sp2可以形成一 σ 鍵，故一個碳 原子會連接 3 個 σ 鍵，而 p 軌域上之電子則會在碳原子所連接的另外三個 碳原子之σ 鍵結上跳動，形成了 π 鍵，這樣的結構稱為 sp2混成軌域，形成 碳-碳雙鍵共價鍵，其楊氏模數與 sp2 _{共價鍵所造成的鍵結強度密切相關，} 所以可以得知原子與原子之間的相互關係對機械性質有一定的影響。 石墨板是由多層石墨所疊成，如圖2.1，層與層之間以 A-B 疊層的方式 組成石墨板結構，如圖2.2。而單層石墨(Graphene sheet)之碳原子以六角形 狀結構規則排列，碳-碳間鍵結距離為 1.42 Å，鍵結角度為 120°。石墨板結 構之單位代表晶體為一六面晶體，相鄰層與層間距(沿著 3 方向即厚度方向) 為3.4 Å。[3]

2.2 原子間之勢能函數

在模擬分子運動時，原子間會有能量，而使原子產生位移、速度、加 速度等各種物理現象，這些物理量都可用勢能函數的形式表示。 在石墨板結構內原子間的作用力的形式可分為五種： (1) 共價鍵結力(Stretch bond)：兩顆原子間距離長度變化而造成之力，如

圖2.3(a)。 (2) 角度鍵結力(Angle bond)：由三顆原子形成一鍵角，當鍵角大小變化 時，會產生一彎矩之內力，如圖2.3(b)。 (3) 扭轉力(Torsion)：由四顆鄰近原子所形成的三個共價鍵，架構成兩平 面所夾之夾角，由角度改變而產生扭轉之力，如圖2.3(c)。 (4) 反轉力(Inversion)：以一顆原子為中心，與另外兩顆原子所組成之平 面，以及第三顆原子所連成之共價鍵，則其平面與共價鍵所夾的夾角 改變而產生反轉力，如圖2.3(d)。 (5) 非鍵結間之力(Non-bonded)：非鍵結力大部分屬於凡得瓦力(van der Waal bond)與靜電力(Electrostatic)，其能量通常隨著任兩顆原子間的距 離增加而迅速減少。 以上均可用勢函數表示，其總能表示如(2.1)式，而本篇論文中選用的勢能 為AMBER force field potential [1]。

bonded non inversion torsion angle bond total U U U U U U = + + + + _{− (2.1)}

2.3 分子力學理論分析

2.3.1 平面載重分析

因為石墨板在受到平面方向之力時，會改變碳與碳間鍵結的距離與角 度，故石墨板的分子力學理論分析當中，推求平面方向(1、2 方向)之材料

參數E1、ν 、12 G12。時，選用勢能為鍵結力勢能以及鍵角勢能之疊加。

(

₀

)

2 2 1 r r k U_stretch = _r − (2.2)

(

)

2 0 2 1_{k θ θ} U_angle = _θ − (2.3) 其中，kr=938 kcal/mol-Å2，kθ=126 kcal/mol-rad2 [23]。 本篇論文是參考Xiao et al.[4]之推導方法，假設在微小變形為下，並簡 化石墨模型，計算模型變形分析以及作用力之影響。 首先，選取代表性的分析單元，並推導碳原子間鍵結長度與鍵結力以 及鍵結角度變化與角度鍵結力間之關係式。將勢能分別對鍵結長度變化量 (△r)與角度變化量(△θ)微分一次得到之方程式即為變形與力之關係式。 r k r F(Δ )= _rΔ (2.4) θ θ = _θΔ Δ k M( ) (2.5) 接著給定簡化之石墨模型一微小變形，可得到作用在碳原子上之力與彎矩。 而平面載重又分為兩部分：(1)單軸方向(1、2 方向)拉伸載重，如圖 2.4 (a)、(2)剪力載重，如圖 2.4 (b)。分別可以得到平面方向的楊氏模數(Young’s modulus)、浦松比(Poisson ratio)與剪力模數(Shear modulus)

2.3.1.1 平面單軸方向拉伸載重：楊氏模數、浦松比 選取代表性的分析單元，並在簡化的石墨模型中，如圖2.5，假設在 Y 方向施以均佈載重，則其合力會作用在OB 及 OA 鍵結上，如圖 2.6。其合 力為 ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = 2 cos 1

α

1

σ

t a f _(2.6) 利用力平衡推導，將 OB 鍵結上集中力

f

與F( bΔ )、M(Δα1)的關係分別列 出，並推求出Δb與Δα₁的關係。利用力平衡推導，則可，所產生的變形關係。 在OB 上鍵結力F( bΔ )與集中力

f

之力平衡式可列為如下： ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = Δ 2 sin ) ( _{b f}

α

1 F _(2.7) OB 鍵結上之彎矩平衡式可列為如下：

[

( ) ( )

]

2 cos 2 1 2 1 α α α _{= Δ − Δ} ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛a _{f M M} (2.8) 其中，定義 2 1

2

α

α

=

−

Δ

Δ

(2.9) 將(2.7)式代入(2.8)式即可得到石墨模型中鍵結力與彎矩關係式。

( )

[

( ) ( )

]

2 tan 2 1 2 1 α α α _{Δ − Δ} ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = Δ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛a _{F b M M} (2.10) 再將(2.4)式、(2.5)式與(2.9)式代入(2.10)式，可得到 OB 鍵結長度變化量與 ∠BOA 角度變化量之關係。 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ Δ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = Δ 2 cot 3 1 1 α α θ k b k a _r (2.11) 由石墨模型中幾何關係變化，可得到1 與 2 方向之應變，如圖 2.7。 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ Δ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + Δ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = 2 sin 2 cos 2 2 sin 1 1 1 1 11

_α

α

α

α

ε

a a b (2.12) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + Δ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − Δ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = 2 cos 1 2 sin 2 2 cos 1 1 1 1 22 _α α α α ε a a b (2.13) 均佈載重應力可由合力與作用截面積關係得到。 ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ Δ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = 2 cos 1 2 sin ) ( 2 cos 1 1 1 1 11 _{α α α} σ ta b F ta f (2.14) 其中

t

為石墨層單層厚度，

a

為鍵結間距，

α

1為鍵結角度。

最後由應力應變關係計算出楊氏模數E及浦松比ν12。 11 11

ε

σ

=

E

_(2.15) 11 22 12

_ε

ε

ν

=

−

_(2.16) 2.3.1.2 平面剪力載重:平面方向剪力模數(Shear modulus) 在簡化的石墨模型中，在周圍施以均佈剪力載重，如圖2.8。則其合力 會作用在OB、OC 及 OA 鍵結上。其等效於集中載重為

f

1、

f

2，則剪力τ 與 等效集中載重之關係式[4]如下。 2 sin 1 1 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ =

α

τ

ta f (2.17) 2 cos 1 1 2 ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = α τ ta f (2.18) 在此同樣利用力平衡推導，將集中力

f

1、

f

2與F( aΔ )、F( bΔ )、M(Δα1)的關 係分別列出，並推求出Δa與Δα₁的關係。利用OB 鍵結上之力平衡關係，可 得到外部剪力與鍵結力之關係式[4]。

( )

a F

( )

b F f f ⎟= Δ =− Δ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ 2 cos 2 sin 1 2 1 1

α

α

(2.19) 其中，由模型受力所產生的幾何變形得知，

b

a

=

−

Δ

Δ

_(2.20) 而OB 鍵結上的彎矩平衡式可列為如下：

( )

2 1 1 1 2 2 cos 2 2 sin 2 M Δα α a f α a f ⎟= ⎠ ⎞ ⎜ ⎝ ⎛ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ (2.21) 將(2.17)式與(2.18)式代入(2.19)式，可得到外部剪力與 OA 鍵結力之關係式。

( )

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

⎟

⎠

⎞

⎜

⎝

⎛

+

Δ

=

2

cos

1

α

1

τ

ta

a

F

(2.22) 再將(2.17)式、(2.18)式與(2.19)式代入(2.21)式，可得到石墨模型鍵結力與彎 矩關係式。

( )

(

2

)

1 1 2 sin 2 cos 1 2 α M α α a F b Δ Δ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ (2.23) 再將(2.4)式與(2.5)式代入(2.23)式，可得到 OA 鍵結長度變化量與∠BOC 角

度變化量之關係[4]。 θ α α α k a k b r ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + Δ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = Δ 2 cos 1 2 sin 2 ₁ 1 2 (2.24) 由模型中幾何關係變化，如圖2.9，可得到 1-2 方向之應變。

(

)

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ Δ − Δ + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + Δ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + Δ = 2 sin 2 2 cos 2 cos 1 2 sin 1 1 1 1 2 12 _α α α α α γ a b a a a a (2.25) 利用應力與應變關係，結合(2.22)式與(2.25)式可得到剪力模數G12。 12 12 12

_γ

τ

=

G

_(2.26)

2.3.2 出平面方向載重分析:出平面方向剪力模數

石墨板的分子力學理論分析當中，假設石墨板為在平面方向為等向性 材料，故在出平面方向之材料係數需另外計算。 因在計算出平面方向材料係數時，考慮層與層間之滑移，所以利用非 鍵結作用力之勢能分析，故選用Lennard-Jones 之勢能公式，來計算石墨層 與相鄰層之相互關係。

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

−

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

=

6 0 12 0

4

ij ij vdW

r

r

r

r

u

U

(2.27) 其中

u

為結合能量( u = 0.0556 kcal/mol )，r0為平衡距離(=3.40 Å) [23]，rij為 原子i與 j_{的間距。此勢能函數以原子間的間距來決定相互關係為吸力或斥} 力，而距離增加時，分子間作用力與勢能會迅速減少(趨近於零)，則在計算 時，故僅考慮rij <rcutoff 內之勢能與作用力，即截斷半徑(cutoff)內之原子，如 圖 2.10。因為由圖 2.11 可看出，凡得瓦勢能與凡得瓦力對於原子間距離的 變化量，可得知在超過約10 Å 之後，凡得瓦力的影響幾乎為零，故本研究 所使用截斷半徑為10 Å。 本篇論文著重於出平面方向之剪力係數G₁₃(=G23)，即在勁度矩陣(2.28) 式中的C44的倒數，相等於柔度矩陣(2.29)式中的(4,4)項。 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = 66 55 44 33 13 13 13 11 12 13 12 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C C C C C C C C C C C C C _(2.28)

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − − − − = 12 13 13 3 1 13 1 13 1 13 1 1 12 1 13 1 12 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 G G G E E E E E E E E E S υ υ υ υ υ υ (2.29) 此處， 1 −

= C

S

(2.30) 則 44 13 C G = (2.31)

本篇論文參考Kelly and Duff [5]之推算演練方法。公式推導方法是將勢 能對1-3 方向應變

γ

13微分兩次而得到C44。 首先定義單顆碳原子對於另一顆原子之勢能

[

6 12 6

]

0 6 0 4 − − − = _{j j} vdW ur r r r U (2.32) 其中r 為原子間距。 j 2 2 2 j j j j

x

y

z

r

=

+

+

(2.33)

接著總和單顆碳原子對於鄰近層截斷半徑內之勢能函數[5]。

[

12 6

]

6 0 6 0

4

ur

r

S

S

U

_vdW

=

−

(2.34) 其中

(

)

∑

∑

− _{= +} − = j n j j n j n r d l S 2 2 2 _(2.35) 2 2 j j j

x

y

l

=

+

_(2.36) 式中

j

為單顆碳原子對應到鄰近層之原子數， d 為石墨層厚度(=3.4 Å)。 因為在每層石墨板中單顆碳原子相對於相鄰層之每一顆碳原子的相對 位置分為兩種平衡位置關係，如圖 2.12，其中一種為碳原子上方為石墨板 六角形結構之頂點(即為碳原子)，此種平衡位置定義為

α

原子，另一為碳原 子上方為石墨板六角形結構之中心，此種平衡位置定義為

β

原子。故計算 時需平均兩種平衡位置之勢能來代表石墨板之碳原子。 單顆原子之能量除以在石墨板中所佔有體積得到勢能密度φ如下[5]：

[

]

4 6 12 6 0 6 0 _{r S S} dA ur ₋ =

φ

_(2.37) 其中 d 為石墨層厚度(=3.4 Å)， A為單顆碳原子在石墨層所佔面積(=2.26 Å2)，而Sn 為兩種不同平衡位置之碳原子Sn值之平均。

2

β α n n n

S

S

S

=

+

(2.38) 接著利用勢能密度相等於應變能密度，

[ ] [ ][ ]

ε

ε

φ

T

C

2

1

=

_(2.39) 其中

[ ] [

]

T 12 23 13 33 22 11

ε

ε

γ

γ

γ

ε

ε

=

_(2.40) 對應變

γ

13微分二次可得到C44， 2 2 2 13 2 13 2 44 j j x x C ∂ ∂ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ≅ ∂ ∂ = φ γ γ φ (2.41) 此處可將應變利用原子位置改變關係，如圖2.13，寫為 0 0 13

d

x

x

_j

−

_j

=

γ

_(2.42) 再利用鏈鎖率得到 0 13 1 d xj = ∂ ∂

γ

(2.43) 代入(2.41)式，可得到[5]

0 0 2 6 2 2 12 2 6 0 6 0 2 0 2 2 2 0 44 4 j j j j x j j x x x j x S x S r dA ur d x d C = = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ∂ ∂ − ∂ ∂ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ≅ φ _(2.44) 其中

(

)

∑

+ + + + + + − = ∂ ∂ j n j j j n j n y x d x n n S n x S 2 ) 4 ( 2 2 2 0 2 2 2 2 ) 2 ( _(2.45) 在此處d 與0 0 j

x

_{均為變形前座標。} 最後由(2.44)式求得出平面方向剪力模數G13之值。

2.4 分子動力學模擬

2.4.1 分子動力學 MD (Molecular dynamics) 模擬

分子動力學模擬的基本原理是以古典牛頓力學為基礎，首先建立一個 系統來模擬微觀世界的物理現象，對於符合古典牛頓力學規律的大量粒子 系統，透過積分方式計算各粒子在空間的運動規律與軌跡，然後利用統計 熱力學方法得出該系統在巨觀的物理特性。 而在建立分子動力學系統時，以建立模擬室(simulation box)為主，模擬 室為一孤立系統(isolated system)，內部原子與能量均不會跑出此系統。在控 制模擬室時，分為各種體系，主要控制體系條件為原子數目、溫度、各方 向壓力、體積。同時必須決定模擬室內在各方向之週期性邊界條件，並選 用Berendsen[27]所發表的數值方法求解。分子動力學在研究 NVT 系綜時， 是將整個環境理想化為熱庫(thermostate)，熱庫的定義為一巨大的可逆熱 源，巨大到任何熱交換都無法改變其溫度，模擬系統熱浴(heat bath)於熱庫 中，藉由彼此間能量的交換，再經過一段時間後，系統的溫度會與熱庫的 溫度相同。Berendsen 溫度修正法，藉由乘上一比例常數來修正原子速度進 而控制系統溫度。而研究 NPT 系綜時，比 NVT 系綜多考慮壓力的影響， 利用與熱浴相同的原理，使模擬系統受到外在壓力的影響，加入修正係數 用來修正模擬室尺寸以及原子位置。本研究在分子動力學利用Berendsen 數 值方法模擬分子運動的過程中，原子間的勢能會使原子受力而產生加速

度，再將加速度對時間積分可以求出原子速度，因熱浴的概念，使原子受 到外部熱庫影響，所以必須修正原子速度，故在 2.46 式中利用χ 當作一比 例常數來控制原子在模擬室內部的速度，使系統的能量可逐漸穩定，其中， T τ 為自訂參數可用來控制比例常數χ (2.47 式)，進而可以控制能量收斂效 果。利用2.46 式求得下一步階 ) 2 1 (t+ Δt 之原子速度，可用2.48 式求得當下時 間(t)的原子速度，並利用2.49 式得到原子之位置。 χ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ Δ + Δ − ← Δ + m t f t t t v t t v ) () 2 1 ( ) 2 1 ( (2.46) 2 / 1 0 1 1 ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − Δ + ← T T t T τ χ (2.47) ← _⎢⎣⎡ − Δ + + Δ)_⎥⎦⎤ 2 1 ( ) 2 1 ( 2 1 ) (t vt t vt t v (2.48) ) 2 1 ( ) (t t r(t) tvt t r +Δ ← +Δ + Δ (2.49) 此數值方法符合Verlet leapfrog 積分方法，式中「←」符號表示右式為對時 間積分可得到左式，為不可逆運算，其中v 為速度，t 為時間， tΔ 為步階時 間，f 為原子所受之力，m 為原子質量，T 為瞬間溫度，T0 為目標溫度，r 為原子位置，「 」符號表示為座標向量。而在模擬室外部受到壓力時，因 內部壓力與外在壓力差而造成模擬室體積大小及形狀的改變，而模擬室各 方面壓力則利用 virail 理論公式計算， 在此 Berendsen 引入壓縮率

(Compressibility)的理論，利用目標壓力 Pext與模擬室當下壓力差，來求得模 擬室的體積變化量，為了計算下一步階時間之模擬室之體積大小及形狀， 將上一步的模擬室體積乘上一比例因子η(2.51 式)，即可得到下一步階的模 擬室體積大小及形狀，而τP為一自訂時間比例常數，用來控制η的大小， 讓模擬室不至於過份變形。 ) 1 ( 1 σ τ η ← -β t P_ext − P Δ (2.50) ( Δ (t+ t) ← Ηt)η Η (2.51) 本研究使用 DL-POLY [28]軟體以 fortran 程式來做分子動力學模 擬，，並利用Material Studio [29]軟體建立以及繪出分子模型，模擬分 析的過程中，步驟上大致要注意下列要點： (1) 勢能函數之選擇 ( Potential function ) 決定分子間作用力之勢能函數後，便可以透過計算求得各種力學性 質。在分子動力學的模擬中，其勢能函數的選取對模擬結果影響極大，選 取不同之勢能函數將得到不同之解答，於是勢能之選擇確實佔有很大的影 響。本章節與第二章相同使用AMBER勢能函數。 (2) 系綜之選擇 ( Ensemble ) 在分子動力學模擬過程中，主要控制條件為原子數目、溫度、各方向 壓力、體積。分子動力學在理論上又分為三種系綜，分別是 NVE、NVT、

NPT。NVE 固定模擬室的原子數量、體積、總能量；NVT 固定模擬室的原 子數量、體積、溫度；NPT 固定模擬室的原子數量、三方向平均壓力、溫 度；NST 固定模擬室的原子數量，各方向壓力、溫度。故在模擬的過程當 中，系綜影響模擬過程中模擬室所代表的物理意義，故系綜的選擇在分子 動力學中相當的重要。

(3) 週期邊界條件（Periodic boundary condition）

利用分子動力學模擬來探討物質的各種現象，為了要精確地描述出真 實物理尺寸下的各種性質，因此我們在取出某一個空間作為代表來探討， 藉由分子動力學模擬在利用週期性的重複排列，而組成真實系統，此方法 稱為週期性邊界條件（periodic boundary condition）。而在本論文中，模擬室 (simulation box)即為模型中的代表性單元(primary cell)，週遭將有重複性投 影(image cell)來表示代表性單元的無窮週期，且模擬過程符合最小映像法則 (Minimum Image Criterion)。

(4) 初始條件設定( Initial condition ) 由於所探討的問題為一動態問題，除了要滿足邊界條件，還必須滿足 初始條件。首先將原子模型排列於模擬室內部，並給定已知的參數，如： 原子位置、速度、壓力和溫度等 (5) 截斷半徑( cutoff ) 選定適當的勢能函數及設定模擬室整體邊界條件與初始條件之後，分

子動力學的模擬即可進行。在模擬的過程中，截斷半徑原理是為了節省電 腦計算時間，取其特定範圍內造成有效影響力之原子做計算，此一範圍即 稱為截斷半徑(cutoff)。即超過截斷半徑外之原子，在模擬過程中不影響解 答結果。分子動力學在模擬過程中，為了配合最小映像法則，DL-POLY 限 制模擬室尺寸邊長須大於截斷半徑。

(6) 系統模擬步數與平衡 ( Number of steps and equilibrium )

分子動力學模擬中必須決定平衡之步階總數，過多的平衡步階數會造 成模擬系統求取時間上的耗費，過少的時間步階可能使系統尚未達平衡狀 態，造成模擬上之誤差。故須確認系統已達收斂平衡穩定狀態，始可觀看 模擬結果。

2.4.2 石墨板模型

本研究利用分子動力學分析方法討論石墨板機械性質。在模擬平面方 向機械性質時，分為三種模型來做模擬比較： (1)單層石墨(Graphene sheet)僅考慮鍵結勢能，不考慮原子間的非鍵結勢能。 在分子動力學模擬(Molecular dynamics)中，分析平面方向機械性質時， 以單層石墨作為其中一種模型。模型(1)僅考慮共價鍵結力與角度鍵結力， 且不加入凡得瓦力之效應，在分子動力學模擬中之模擬室為週期性邊界條 件，其 1、2 方向鍵結為連續性的週期性邊界，單層石墨板模擬室尺寸為 47.3Å×52.1Å×3.4Å，由 924 顆碳原子組成。石墨層單層厚度為 3.4Å。如圖

2.14。 (2)單層石墨(Graphene sheet)考慮鍵結勢能以及非鍵結勢能。 在分子動力學模擬(Molecular dynamics)中，模型(2)考慮共價鍵結力與 角度鍵結力，並加入凡得瓦力之效應，在分子動力學模擬中之模擬室為週 期性邊界條件，其1、2 方向鍵結為連續性的週期性邊界，僅考慮單層石墨 的凡得瓦力效應，且截斷半徑為 10Å，故將模擬室厚度設定為 34Å。而石 墨板模擬室尺寸為 47.3Å×52.1Å×34Å，由 924 顆碳原子組成。石墨層單層 厚度為3.4Å，如圖 2.15。 (3)多層石墨(Graphene flake)考慮鍵結勢能以及非鍵結勢能。 在分子動力學模擬(Molecular dynamics)中，多層石墨以八層石墨作為模 型，層與層間之間距為3.4 Å，故將模擬室厚度設定為 27.2Å，考慮共價鍵 結力與角度鍵結力，並加入凡得瓦力之效應。在分子動力學模擬中之模擬 室為週期性邊界條件，其1、2 方向鍵結為連續性的週期性邊界，考慮石墨 層與層間的凡得瓦力效應，且截斷半徑為 10Å。此石墨板模擬室尺寸為 47.3Å×52.1Å×27.2Å，由 7392 顆碳原子組成。石墨層單層厚度為 3.4Å，如 圖2.16。 在石墨模型上施單軸應力或應變來計算平面方向楊氏模數、浦松比。平面 方向剪力模數方面，對石墨模型施以剪應力，經過模擬計算後求得剪應變，

再利用應力-應變關係得到平面方向剪力模數。並利用模型(1)在平面方向的 模擬結果與解析解比較。 本研究利用上述三種模型來得到平面材料參數。而模擬出平面方向機 械性質時，則利用多層石墨模型，即模型(3)，在考慮鍵結勢能與非鍵結勢 能的狀況下，施以出平面方向剪應力來得到出平面方向之剪力模數，此模 擬結果將與解析寫比較。 利用分子動力學模擬石墨板的過程中，分為兩個步驟：首先將模擬室設定 目標壓力0 kbar、目標溫度 0K，利用 NST 系綜進行平衡，模擬步數為 100000 步，每一步階時間為0.001ps，總模擬平衡的時間為 100ps，給定數值方法[27] 中自定參數τ_T =1、τ_P =10。 確定石墨結構產生新的平衡狀態之後，圖 2.17、圖 2.18 分別為模型(1)、 模型(2)在模擬過程中(a)能量變化圖、(b)溫度變化圖、(c)1 方向壓力變化圖、 (d)2 方向壓力變化圖、(e)3 方向壓力變化圖；圖 2.19 為平衡穩定狀態下， 壓力為零之石墨板。接著將上一步驟所產生新的模擬室狀態當作初始狀 態，同樣利用NST 系綜，對模擬室施以應力或應變模擬而達到平衡收斂， 產生新的石墨受力平衡狀態。而在模擬的過程中，總能量、溫度及壓力各 項將會快速地達到收斂，主要原因有三個，其中之ㄧ為石墨本身即為一穩 定簡單之結構，故在平衡時不易產生原子大幅度跳動；其二是因為分子動 力學計算過程中的設定前5000 步，為防止原子過度移動而強制將分子速度

乘上一比率，使原子慢慢地趨於穩定而使整體系統平衡；其三為數值方法 中自定參數的改變會影響收斂性之快慢，而選用適當的參數則會有良好的 收斂效果。

2.4.3 石墨板平面方向機械性質

2.4.3.1 平面單軸方向拉伸載重：楊氏模數、浦松比 此章節計算楊氏模數與浦松比時，分別對石墨模型施予應力及應變兩 種方法來求得楊氏模數及浦松比，並比較兩種方法所得到之結果是否接近。 (1)施予石墨板單方軸應力(axial stress)得到楊氏模數。 石墨板模型的LB、WB分別為模擬室(simulation box)長度及寬度、h 為 厚度，如圖2.14、圖 2.15、圖 2.16。設定目標溫度為 0K 之狀態，利用 NST 系綜進行平衡，同樣地，模擬步數為100000 步，每一步階時間為 0.001ps， 模擬時間為100ps，給定數值方法中自定參數τ_T =1、τ_P =20，在拉伸方向(1 方向)對模擬室施以應力 σBox = 0.1 kbar，如圖 2.20，此模擬條件適用於模型 (1)~ (3)，最後確認新的模擬室狀態達到平衡收斂，如圖 2.21。由於模型(2) 的 石墨板實際所得之應力與模擬室上不同，故須經過比例換算，以得到石墨 所受實際應力。 number layer × = Graphene Box Box Graphene t h σ σ _(2.52)

其中，σGraphene為石墨所受實際應力，σBox為模擬室所受外部應力，而tGraphene 為單層石墨厚度(3.4Å)。經過 MD 對內部原子平衡後，可得到模擬室在拉伸 方向及橫向方向(2 方向)之應變。 B B L L Δ = 11 ε _(2.53) B B W W Δ = 22 ε (2.54) 利用應力-應變關係，即可求得石墨板之楊氏模數。 11 2 1 _ε σB E E = = (2.55) 另外，模擬室橫向方向之應變，與拉伸方向之應變之比值，即為石墨模型 的浦松比。 11 22 12 _ε ε υ =− (2.56) (2)施予石墨板均勻單方向應變，利用勁度矩陣求得楊氏模數。 首先，設定目標溫度為0K 之狀態，利用 NVT 系綜進行平衡，同樣地， 模擬步數為 100000 步，每一步階時間為 0.001ps，總模擬平衡的時間為 100ps，給定數值方法中自定參數τ_T =1、τ_P =20，在拉伸方向(1 方向)對模 擬室施以應變0.001，如圖 2.22，此模擬條件用於模型(2)~ (3)，最後確認新 的模擬室狀態達到平衡收斂，並得到模擬室各方向應力。由於模型(2)的石

墨板實際所得之應力與模擬室上不同，同樣須利用 (2.52)式經過比例換 算，以得到石墨所受實際應力。透過平面應力-應變關係可得勁度矩陣中 11 C 、C₁₂、C₂₁、C₂₂之值，如表2.1，並推導出石墨板之楊氏模數及浦松比。 ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ 12 22 11 66 22 21 12 11 12 22 11 0 0 0 0

γ

ε

ε

τ

σ

σ

C C C C C (2.57) 22 21 12 22 11 C C C C C E = − (2.58) 22 12 12 C C = ν (2.59) 2.4.3.2 平面剪力載重:平面方向剪力模數(Shear modulus) 首先，設定目標溫度為 0K 之狀態，利用 NST 系綜進行平衡，模擬步 數為200000 步，每一步階時間為 0.001ps，總模擬平衡的時間為 200ps，給 定數值方法中自定參數τT =1、τP =20，對模擬室施以 1-2 方向剪應力 τBox = 0.001kbar，使石墨板產生 1-2 方向之剪應變，如圖 2.23、圖 2.24 為模型(1)~ (3)施以 1-2 方向之剪力示意圖。此模擬條件適用於模型(1)~ (3)，由模型(2) 石墨板所受之實際應力與模擬室上不同，故須經過比例換算，以得到石墨 板所受實際應力。

number layer × = Graphene Box Box Graphene t h τ τ _(2.60) 其中，τGraphene 為石墨所受實際剪應力，τBox為模擬室所受外部應力，經過 能量平衡之後，模擬室產生1-2 方向之剪變形。將石墨受到實際剪應力與剪 應變相除，可求得剪力模數。 12 12 _γ τ_Graphene G = _(2.61)

2.4.4 出平面(Out-of-plane)方向載重分析:出平面方向剪力模數

此模擬條件適用於模型(3)，首先，設定目標溫度為 0K 之狀態，利用 NST 系綜進行平衡，模擬步數為 500000 步，每一步階時間為 0.001ps，總 模擬平衡的時間為 500ps，給定數值方法中自定參數τT =1、τP =20，則對 模擬室施以1-3 方向 0.001kbar 之剪應力，如圖 2.25，此模擬條件適用於模 型(3)，最後經過能量平衡之後，模擬室產生 1-3 方向之剪變形，利用應力 應變關係求得剪力模數。 13 13 13 _γ τ = G _(2.62)

2.5 結果與討論

2.5.1 石墨板材料模數與參考文獻之比較與討論

此研究利用解析解與分子動力學兩種方法研究石墨板的材料模數，碳 與碳之間之作用力引用 AMBER[1]的鍵結勢能之公式參數和 T.- W. Chou [23]所採用的 Lennard-Jones 非鍵結勢能公式參數來分析石墨板機械性質。 而石墨板分為解析解與分子動力學解兩種分析方法，並討論解析解與分子 動力學解兩者解答的接近性。 表2.1 中各項勁度係數表示此為一對稱矩陣，可顯示出石墨板在平面方 向為等向性材料。表 2.2、表 2.3 分別為平面方向楊氏模數、浦松比和平面 方向剪力模數之解析解與分子動力學解分析之結果，用解析解求得之楊氏 模數為0.805 TPa，另外，分子動力學模擬：模型(1)利用應力-應變法求得楊 氏模數的解答為0.794 TPa。結果顯示在平面載重分析方面，解析解與模型 (1)的解答誤差約在 1.3 %以內，所以相當的接近。在楊氏模數方面，利用應 力法以及勁度矩陣法所得之楊氏模數及浦松比結果相同，而利用分子動力 學模擬考慮凡得瓦力之石墨板均比不考慮凡得瓦力的石墨板楊氏模數要高 出一些，其中，模型(2)的楊氏模數又比模型(3)高出許多，表示只考慮單層 石墨板原子間的凡得瓦力，結果會有明顯的偏高。 參考文獻所模擬的結果介於0.671~1.153 TPa 之間，本研究的楊氏模數 結果均介於此範圍內，且文獻[7]中顯示使用 AMBER 勢能分析石墨板平面

方向材料模數時，會比使用 Morse 勢能時低，故使用 AMBER 勢能所模擬 之楊氏模數在合理的範圍之內。而本篇計算浦松比結果為0.273，其他文獻 的浦松比分析結果介於0.195~0.428 之間，所有理論分析的值均高於實驗值 [9]。 而 平 面 方 向 剪 力 模 數 的 計 算 結 果 ， 參 考 文 獻 的 剪 力 模 數 介 在 0.384~0.482 TPa，本研究所模擬之結果稍微低於參考文獻的值。 本研究利用分子動力學以及解析解兩種方法，分別探討之楊氏模數、 浦松比和剪力模數，其結果均符合等向性材料G12 = E1/2(1+ν12)之特性。 表2.4 為出平面方向剪力模數之解析解與分子動力學解分析之結果，用 解析解得到出平面剪力模數為0.209 GPa，而用分子動力學模擬得到出平面 剪力模數為0.290 GPa，因為在計算凡得瓦勢能所考慮的石墨層數不同，因 此會有較大的差異。其他文獻結果的範圍在0.18~0.42 GPa 之間，表示理論 結果與勢能之選用有很大的關聯。 本研究中，利用解析解與分子動力學模擬兩種方法所得到之石墨板平 面方向材料參數誤差在0.6~2.9%之間，而出平面方向材料參數誤差 38.7%， 顯然用分子力學理論分析與分子動力學模擬兩種分析方法來預測石墨板的 機械性質，可以得到相當接近的結果。 在平面方向與出平面方向的材料參數分析中，其相鄰層石墨均以凡得 瓦力(van der Waals bond)互相吸引連結，因為單層石墨層中碳原子相互間之

鍵結勢能比石墨層與層之間的凡得瓦力要來的大，故此現象驗證在材料參 數中，可以發現平面碳板強度相較出碳板強度高出許多。 由於出平面方向剪力模數較小，則碳原子在受力之後的相對位置上， 石墨層與層之間的移動與平面方向碳原子之角度或間距改變相較之下，更 容易產生石墨層之間的滑動。所以石墨板在受力時，石墨層與層之間的滑 動會比平面方向的拉伸或剪力變形更容易產生。 本研究分析結果與實驗結果比較，在平面方向材料模數有相當好的比 對結果，而在出平面方向材料模數則有較大的誤差。在目前不同的研究文 獻中，出平面方向之剪力模數方面，其值約在 0.13~4.5 GPa 相當大的範圍 內，因為沒有較固定之值，所以無從判斷其正確之理論值，僅與實驗值時 做比較。

2.5.2 不同凡得瓦力勢能對於模擬之結果及影響

在石墨板的研究當中，最後又引用 P.M. Agrawal、I.M. Daniel、Dreiding 所發表的凡得瓦力勢能公式及參數，整理於表2.5，並將各種不同的勢能公 式參數轉換為同樣的形式及單位。由表中可得知各種勢能有不同的碳-碳平 衡間距r₀和結合能量，可預估出模擬結果將會有不同的數據結果。 如表 2.6~2.8，同樣利用模型(2)~(3)中，相同尺寸下之模擬室，使用分 子動力模擬方法計算出石墨板的各種材料模數，而模擬過程中石墨板的截 斷半徑均為10Ǻ。