修士論文

分子動力学法によるカーボンナノコイルの構造・力学特性評価

指導教員： 田中

毛利

克志

友宙

2015 年 2 月
神戸大学大学院 工学研究科

博士課程前期課程

機械工学専攻Master thesis

Molecular Dynamics Studies
on Structure and Mechanical Characteristic ofCarbon Nano CoilTomohiro MOHRI
Feburary 2015
Department of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, Japan要約
CNC はナノスプリングや電磁波吸収材としての応用が期待される材料である．本研
究では，単層及び複層の CNC に対して分子動力学シミュレーションを用いて初期構
造緩和，引張，繰り返し変形のシミュレーションを行うことで CNC の機械的特性について検討した．Zigzag，Armchair，Chiral 型の 3 種類のパターンの短いストレート
チューブ切片を 36 角形状に接合することで CNC を作成しており，熱処理を与えるものと与えないものの二つの条件で構造緩和を行うことで，構造安定性について検討した．コイル径が 5[nm] の小規模な CNC は，熱処理を与えない条件ではどのモデルもら
せん構造が安定しなかった．熱処理を与える条件では Chiral 型のみらせん構造が安定した．コイル径が 20〜30[nm] の大規模な SWCNC に対しては，熱処理を与える条件についてのみ構造緩和を行い，安定した CNC に対して引張シミュレーションを行った．
いずれの CNC も，λ=1.5 程度までは線形的に応力上昇しており，ばね定数は 10−3 程度のオーダーを示していた．一部，応力が二次曲線的に上昇する硬化挙動を示す CNCも存在した．どの SWCNC もチューブに局所的な座屈やねじれが生じることで応力が減少した．安定した SWCNC に対して引張後除荷を行う繰り返し変形シミュレーショ
ンを行った結果，全ての SWCNC で除荷開始時に負荷反転時の慣性の影響により応力に大きな揺らぎが生じたが，除荷が進むと揺らぎは減少し除荷時の応力-ストレッチ曲
線は負荷時と同じ経路に戻った．多くの SECNC で負荷時に発生した局所的な座屈が除荷時に解消されており，チューブ変形の復元性が高いことが確認出来た．安定した

SWCNC を重ねることで 2 層の MWCNC を作成し，構造緩和シミュレーションを行った結果，グラファイトの層間隔である δ=0.35nm よりも層間距離が大きい MWCNC では安定したらせん構造は得られなかったが，層間距離が δ=0.35nm であるカイラリティ(30，20)-(40，20) と (40，20)-(45，25) の MWCNC はらせん構造が安定した．らせん構造が安定した MWCNC に引張シミュレーションを行った結果，いずれもばね定数は2.65×10−3 [µN/nm] 程度を示し，ピーク応力後の応力減少はチューブにねじれによって局所的な座屈を生じることによってもたらされていた．

Summary
Carbon nano coil(CNC) is expected as nanospring and electromagnetic wave absorbent. In this study, various molecular dynamics (MD) simulations are performed on
the single walled CNCs(SWCNCs) and multi walled CNCs(MWCNCs). We made various CNCs and investigated their stability with/without heat treatment, The CNCs are
made by connecting 36 straight segments of zigzag, armchair and various chiral carbon
nanotubes(CNT) in the coil radious of 5nm. The results show that no CNC keep thecoil form without heat treatment, and only the chiral CNCs keep the coil form after
the heat treatment. Then we made larger chiral CNCs of the coil radius 20〜30nmwith the heat treatment, and performed tensile simulations on the stable CNCs. Most
CNCs show linear stress-stretch response up to the stretch of λ=1.5, and the spring
constant is in the order of 10−3 [µN/nm]. Some CNC shows the strain hardening, resulting in the quadratic curve of the stress-stretch relation. All CNCs show stress dropby the local buckling and twist of the tubes. Then we performed cyclic loading on thechiral CNCs. Due to the inertial motion in the inverse of the loading, almost all the
CNCs show the large stress oscillation in the early stage of the unloading; however,
the oscillation diminishes and ﬁnally the unloading stress-stretch relation converges tothat in the loading process. Observations on the CNC deformation reveal that many
CNCs show recovery of the local buckling in the cyclic loading. Finally we performed
simulations about the stability and tensile behavior of nested two walls CNCs. The
nested CNCs of the segment chirality (30,20)-(40,20) and (40,20)-(45,25) are stable,
and nested spacing is same as the stable distance of the graphite layer, δ=0.35nm. On
the other hand, the other two tubes with the larger tube spacing couldn t keep thecoil form. The tensile simulations on these stable MWCNCs revealed that the spring
constant is about 2.65×10−3 [µN/nm] for both, and the stress-peak is brought by the
local buckling and twist.

目次
第1章 緒 論1

第 2 章 解析手法の基礎

52.1

分子動力学法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52.2

原子間ポテンシャル. . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.2.1

Brenner ポテンシャル . . . . . . . . . . . . . . . . . . . . . . .

6

2.2.2

Lennard-Jones ポテンシャル . . . . . . . . . . . . . . . . . . . .

10

2.2.3

力の表式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.3

高速化手法

2.4

速度スケーリング法. . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.5

カーボンナノチューブの幾何形状パラメータ . . . . . . . . . . . . . . .

17. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 3 章 単層 CNC の構造安定性

3.1

3.2

14

19

小規模解析モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

3.1.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

3.1.2

解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

大規模解析モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.2.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.2.2

解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

第 4 章 単層 CNC の引張シミュレーション

33

4.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

4.2

解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

第 5 章 単層 CNC の繰り返し変形シミュレーション

58

5.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

58

5.2

解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

58

i

目次
第 6 章 2 層 CNC の構造安定性・引張シミュレーション

6.1

6.2

ii
75

初期構造緩和シミュレーション . . . . . . . . . . . . . . . . . . . . . .

75

6.1.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

6.1.2

解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

76

引張シミュレーション . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

6.2.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

6.2.2

解析結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

第7章 結 論

86

参考文献

89

学術論文・学術講演

91

謝辞

94

第1章
緒論
カーボン材料として黒鉛やダイヤモンドが古くから知られていたが，近年，Kroto，

Smalley，Curl らによって発見されたフラーレン [1] や，飯島らによって発見されたカー
ボンナノチューブ (Carbon nano tube : CNT)[2] など，ナノメートルレベルでの新たな材料が次々と発見され，これらの物性についての研究が盛んに行われている．これ
らのナノカーボン材料の一つにカーボンナノコイル (Carbon nano coil : CNC) と呼ばれる材料がある（図 1.1）[3] ．CNC は CNT をばねのようならせん状に巻いたような構造を持っており，コイル外径が数十から数百 [nm] 程度の炭素繊維である．CNC の特徴としては，(1) 極細のらせん状炭素繊維である，(2) 極小のコイルばね構造である，(3) 電磁波吸収特性に優れる，(4) 低い電界で電界電子放出する，(5) 伸縮性がありばねのように振る舞う [4] といったものがある．これらの優れた特性から電磁波吸収材としての応用が期待されるだけではなく，ナノ電気機械システム分野におけるナノイン
ダクタやナノスプリングとしての応用などを期待した研究がおこなわれている．また，

CNT はその直径とグラフェンシートのねじれ具合から，性質が半導体的であったり金属的であったりすること [5] が知られているが，CNC では，二つに加えて半金属状態の性質を示す様態の存在が指摘されており，超電導コイルの可能性が期待されている [6] ．

1950 年代に Davis，Slawson，Rigby らが CNC の合成 [7] について初めて報告を行っており，ナノカーボン材料の中では比較的初期に登場したが，偶発的に発見されたものであり有効な合成条件は見つかっていなかった．そのため，CNC を具体的な対象とした研究は盛んではなかったが，1975 年に Baker らが再現性の高い CNC の合成手法を確立し [8] ，1990 年に元島らがカーボンマイクロコイル (Carbon micro coil : CMC) の

1

第1章

緒論

2

合成を報告する [9] など徐々に CNC に対する関心が高まってきた．そして，1991 年に

CNT が発見され，その高い特性から CNT の研究が非常に盛んに行われるようになると同時に，直線的な CNT だけでなく，らせん構造を有した CNC についても高い関心が集まるようになった．
これまで，カーボンナノ材料に対する研究は多く報告されており，特に CNT については非常に多くの研究が行われている．実験によるアプローチでは，CNT を発見した
飯島らがさらに 1993 年に単層 CNT(Single walled Carbon nano tube : SWCNT) の合成に成功しており [10] ，Nasibulin らが触媒微粒子の粒径と SWCNT の直径の相関について報告している [11] ．Kawakubo らは，複層の CNT(Multi walled Carbon nano tube

: MWCNT) の粉末を撒いた平面摩擦実験を行い，MWCNT 粉末のチューブ径が大きいほど摩擦係数が小さいことを示した [12] ．阿知波らは，レーザーアブレーション法における実験条件を精密に調節することで，特定のカイラリティからなる SWCNT の合成に成功した [13] ．小山らは燃料電池の触媒担体に CNT を用いることで燃料電池の性能が向上すると報告している [14] ．計算材料科学の分野でも盛んに研究は行われており，
尾方，渋谷らは五員環や七員環のような欠陥 (Stone Wales 欠陥) を持った CNT について，熱伝導率の評価を分子動力学法を用いて行っている [15] ．渋田，丸山は Fe,Co,Niの 3 つの触媒金属において，触媒金属種の違いが CNT の生成過程に与える影響について分子動力学法を用いて検討している [16] ．Hirai らは，ピンホール欠陥を持つ CNTの機械的特性について調べており [17] ，Deguchi，Yamaguchi らは，SWCNT に高温下で引張応力を与えることで Stone-Wales 欠陥を発生させ，引張時におけるそうした欠陥の挙動について評価している [18] ．西村らは，MWCNT の解析として，2 層構造を持つ CNT の圧縮特性について分子動力学法によって検討している [19] ．

CNC に関する実験研究として，Pan，Zhang，Nakayama らは熱 CVD 法による CNCの安定合成法 [20] ，細川，後藤らは CVD 法による CNC の大量合成法を報告している[21]．Nakayama，Pan らは，アセチレンの流量や触媒を調整することで CNC のコイル

径，長さの制御に成功している [22] ．Katsumata らは CNC と CNC よりも捻じれの強いカーボンナノツイスト (Carbon nano twist : CNTw) の作り分けが可能であると報告し [23] ，久米，長谷川らは，CNC の先端触媒粒子を TEM 観察することで触媒粒子からのカーボンの排出速度の差が CNC の成長機構の一つであると報告している [24] ．計

第1章

緒論

3

算材料科学の分野でも，香川，山口，平原らが分子動力学法を用いて CVD 法による

CNC の生成過程の再現を試みており，CNC の成長機構の解明に取り組んでいる [25] ．
Lizhao，Jijun らは CNT を六角形状に繋ぐことで CNC のモデルを作成し，CNC のコ
ンダクタンスについて検討している [26] ．他にも CNC に関する研究は多く存在するが，CNT と比較するとまだまだ少なく，特に計算分野では CNC の詳細な構造がまだ解明されておらず，モデルの作成が難しいため報告例が少ないのが現実である．また，作成されたモデルも実際に存在する CNC と比較すると非常に規模の小さなものがほとんどである．
本研究では，様々なカイラリティのストレートチューブ切片を 36 角形状に接合することで CNC のモデルを作成し，分子動力学シミュレーションにより初期構造緩和，引張，繰り返し負荷などを行うことで CNC の構造安定性ならびに力学特性について解明することを主たる目的とする．また，作成した SWCNC を組み合わせることで 2 層の MWCNC も作成し，検討を行う．
第 2 章では解析手法の基礎として，分子動力学法を簡単に説明し，分子動力学計算で最も重要となるポテンシャルエネルギーについて述べる．また，大規模計算を行うための高速化手法を示し, 最後に CNT 特有のカイラリティについて説明する．
第 3 章では，コイル径やカイラリティの異なる様々な単層 CNC に対して，熱処理を与えるものと与えないものの二つの条件で構造緩和を行い，CNC の幾何形状の変化か
ら構造安定性について検討を行う．
第 4 章では, 第 3 章で安定した SWCNC を対象に引張シミュレーションを行う．2 通りのひずみ速度で引張を行い，SWCNC の幾何形状やひずみ速度が力学応答やばね定数に与える影響について議論する．
第 5 章では，第 4 章の引張の後に負荷ひずみ増分を反転させて元の周期長さまで戻す
繰り返し変形シミュレーションを行い，負荷時と除荷時の応力挙動の違いや SWCNCの復元性について議論する．
第 6 章では，第 3 章で安定した SWCNC のうちから数種類の組み合わせで 2 層のMWCNC を作成し，構造緩和シミュレーション及び引張シミュレーションを行うこと
で，単層と複層での幾何形状やばね定数，応力応答の違いについて議論する．
最後に, 第 7 章で本研究の総括を述べる.

第1章

Fig.1.1 SEM image of CNC[3]

緒論

4

第2章
解析手法の基礎2.1

分子動力学法

分子動力学法 (molecular dynamics method，略して MD 法) は，系を構成する各粒子についてニュートンの運動方程式

mid2 r i
= Fi
dt2

(2.1)

を作成し，これを数値積分することにより粒子の軌跡を求める方法である．ここで，mi ，r i はそれぞれ粒子 i の質量および位置ベクトルである．粒子 i に作用する力 F i は，系のポテンシャルエネルギー Φtot の各位置における空間勾配として次式により求められる．

Fi = −

∂Φtot
∂r i

(2.2)

式 (2.1) の数値積分には，Verlet の方法，予測子–修正子法等がよく用いられる．本
研究では，以下に示す Verlet の方法を用いた．
時刻 t + ∆t と t − ∆t での粒子 i の位置ベクトル r i (t ± ∆t) を Taylor 展開すると，

dr i (t) (∆t)2 d2 r i (t) (∆t)3 d3 r i (t)r i (t + ∆t) = r i (t) + ∆t
+
+
+・
・
・
dt
2
dt2
3!
dt3
dr i (t) (∆t)2 d2 r i (t) (∆t)3 d3 r i (t)
+
−
+・
・
・
r i (t − ∆t) = r i (t) − ∆t
dt
2
dt2
3!
dt3

(2.3)
(2.4)

となる．ここで，v i を時刻 t における粒子 i の速度とすると，

dr i= v i (t)
dt
5

(2.5)

第2章

解析手法の基礎

6であり，式 (2.1) と式 (2.5) を式 (2.3) と式 (2.4) に代入すると，(∆t)2 F i (t) (∆t)3 d3 r i (t)r i (t + ∆t) = r i (t) + ∆tv i (t) +
+
+・
・
・
2
mi
3!
dt3
(∆t)2 F i (t) (∆t)3 d3 r i (t)r i (t − ∆t) = r i (t) − ∆tv i (t) +
−
+・
・
・
2
mi
3!
dt3

(2.6)
(2.7)

となる．両式の和と差をとると，

F i (t)+・・・mi(∆t)3 d3 r i (t)+・・・r i (t + ∆t) − r i (t − ∆t) = 2∆tv i (t) + 2
3!
dt3
r i (t + ∆t) + r i (t − ∆t) = 2r i (t) + (∆t)2

(2.8)
(2.9)

が得られる．これより，時刻 t + ∆t での位置ベクトルと t での速度は

r i (t + ∆t) = 2r i (t) − r i (t − ∆t) + (∆t)2v i (t) =F i (t)mi1{r i (t + ∆t) − r i (t − ∆t)}2∆t(2.10)(2.11)と求められる．t + ∆t での座標を求めるには 2 つの時刻 t と t − ∆t での座標が必要である．初期の計算 (t = 0) では，t = ∆t での座標 r i (∆t) は式 (2.6) と初速度から得ることができる．

2.2

原子間ポテンシャル

粒子に作用する力は系のポテンシャルエネルギ−により決定される．そのため，分子動力学法においてはポテンシャルの選定が重要になる．解析の対象となる物質や解析条件に適当となるようなポテンシャル関数を決定しなければならない．

2.2.1

Brenner ポテンシャル

本解析では，シリコンに対する Tersoﬀ 型ポテンシャルが Brenner によりグラファイ
トにフィッティングされたものを用いる [27] ．
シリコンは常温常圧においてダイヤモンド構造を持つが，炭素の場合ダイヤモンドとグラファイトという 2 つの安定構造がある．したがって，グラファイトの sp2 結合と，ダイヤモンド構造の sp3 結合の違いを表現することが重要になる．

第2章

解析手法の基礎

7

系のエネルギーはΦtot =∑ ∑[]fc (rij ) VR (rij ) − B¯ij VA (rij )(2.12)i j(>i)と表される. rij は粒子 i, j 間の距離を示している．ここで VR (rij ) は斥力をあらわす項であり，−B¯ij VA (rij ) は引力をあらわす項である．また fc (rij ) はポテンシャルエネ
ルギーの打ち切り（カットオフ）を滑らかにするための項である．各項はそれぞれ[]√De
exp −β 2S (rij − Re )
S−1


√
S · De2
VA (rij ) =exp −β(rij − Re )
S−1
S

VR (rij ) =



1,

 [


fc (rij ) =






(

1 + cos

π (rij − R1 )R2 − R1(2.13)(2.14)
rij < R1)]

/2, R1 < rij < R20,

(2.15)

rij > R2である. B¯ij は，粒子 i ，粒子 j 以外に粒子 k も含めた 3 個の原子によって定められる項であり,

1B¯ij = (Bij + Bji )2(2.16)Bij , Bji はそれぞれ原子 i, j を中心とした加算となり一般に Bij ̸= Bji である.
Bij = [1 +∑G(θi )fc (rik )]−δ(2.17)

G(θj )fc (rjk )]−δ(2.18)

k̸=i,jBji = [1 +∑k̸=i,j

ここで rik は粒子 i, k 間の距離, rjk は粒子 j, k 間の距離を示している. θi は原子 i を中心とする j − i − k の内角, θj は原子 j を中心とする i − j − k の内角である. 微分のため,∑G(θi )fc (rik ) を ζij ,∑G(θj )fc (rjk ) を ζji とする. G(θ) は
[

c20
c20
G(θ) = a0 1 + 2 − 2
d0 d0 + (h + cos θ)2](2.19)

である.
式 (2.12) ∼ 式 (2.19) で用いられるパラメータは C2 分子，独立したグラファイト
シートおよびダイヤモンドそれぞれの結合エネルギと平衡状態における結合距離，仮第2章

解析手法の基礎8

想的な単純格子および f cc 構造の結合エネルギ，グラファイトの斜方六面体構造から
ダイヤモンド構造への相変態における障壁エネルギの計算等によりフィッティングされている．それらのパラメータを表 2.1 に示す．

Table 2.1 Potential parameters for Brenner potential.
De [eV]
6.0
Re [nm]0.139
β [nm−1 ]
0.21
S
1.22
h
1.0
a0
0.00020813
c0
330.0
d0
3.5
R1 [nm]0.17
R2 [nm]0.20

シリコンに対するパラメータの場合と最も異なることは，h の値が 1 になっている点である．h が 0 の場合はより小さな結合角で安定するのに対し，h を 1 にすることで結合角はより大きく，すなわち 180 [deg.] に近づこうとする．このことが，グラ
ファイトシート内の平面構造を維持する駆動力になる．また Brenner ポテンシャルには, 炭素原子間距離の値に重点を置きクラスタの形成に最適化されたパラメータ 1 と,炭素間に作用する力の値に重点を置き物性の測定に最適化されたパラメータ 2 が存在
する. 本研究では力の再現を重視したパラメータ 2 を用いて計算を行った.
結合距離の変化に対する各結合角でのポテンシャルおよびポテンシャルが最小の値をとる時の結合角の変化に対する結合距離の変化を図 2.1，2.2 に示す．図 2.1 より, ポ
テンシャルの最小値が結合角の減少に伴い結合距離が大きくなる方向へ推移しているのが確認できる. また図 2.2 より結合角 60∼90[deg.] において, 結合距離が増減し変化が滑らかではないのは, 原子 j, k での結合長に対するカットオフ半径による.

第2章

解析手法の基礎

Potential energy E [eV]20

θ = 30θ = 60θ = 120
θ = 180
10

0

0.1

0.12

0.14
0.16
0.18
Bond length r [nm]0.2

Fig.2.1 Relationship between potential energy and bond length.

Stable bond lengthr [nm]

0.2

0.18

0.16

0.14
60

120

180

Bending angle θ [deg]Fig.2.2 Relationship between stable bond length and bending angle.

9

第2章

2.2.2

解析手法の基礎

10

Lennard-Jones ポテンシャル

グラファイトの層間の Van der Waals ポテンシャルは，次式のような Lennard-Jones型ポテンシャルで表される [28] ．
∑ ∑

(

σ
Φtot =4ϵ riji j(̸=i)

)12

(

σ
−
rij

)6 


(2.20)

第 1 項は斥力，第 2 項は引力を表し，式 (2.20) で用いられるパラメータを表 2.2, ポテ
ンシャル曲線を図 2.3 に示す．

Table 2.2 Potential parameters for Van der Waals.

ϵ
σ0.004783 [eV]0.3345 [nm]–3

Van der Waals potential E [eV]

×10

5.0

0.0

–5.0
0.3

0.4
Atomic distance

0.5

0.6

r [nm]

Fig.2.3 Relationshipbetween Van der Waals potential and atomic distance.
r

第2章

2.2.3

解析手法の基礎

11

力の表式

本解析で用いる各ポテンシャルについて，式 (2.2) に基づいて各ポテンシャル関数
より原子間力を導く．まず，Brenner ポテンシャルについて原子 ij 間のエネルギー寄与は

Φij = VR (rij ) −Bij + Bji
VA (rij )2(2.21)

である. これの r i , r j , r k による微分

Fi = −

∂Φij
∂Φij
∂Φij
,Fj = −
,Fk = −
∂r i
∂r j
∂r k

を, j > i 全ての組について加算すれば原子に働く力が求まる. 以下で個々の原子にかかる力を示す.
原子 i の力は]r[[∂Bij ∂Bji1
F i = B¯ij VA′ (rij ) − VR′ (rij )+ VA (rij )
+
rij
2
∂r i
∂ri

]

ij

(2.22)

である. ここで ′ は r による微分を表す. VR (rij ), VA (rij ) の r ij での微分は
]r
[
√
√
De
ij
VR′ (rij ) = −β 2S
exp −β 2S (rij − Re )S−1
rij

(2.23)

2 S · De2r ij
exp −β(rij − Re )
S S−1
S
rij

(2.24)√VA′ (rij ) = −β





√

である. fc (rij ) の r ij での微分は
(



 −1

∂fc (rij )
=

∂r ij


π
π (rij − R1 )sin
2 R2 − R1
R2 − R1)0: R2 > rij > R1

(2.25): rij < R1 , R2 < rijである. 以降では Bij , Bji の r i での微分を考える.

∂Bij
∂ζij
= −δ(1 + ζij )−δ−1
∂r i
∂r i

(2.26)(

中心∂ζij i∑r ik∂ cos θi=fc′ (rik )G(θi )+ fc (rik )G′ (θi )
∂r i
rik
∂r i
k̸=i,j)(2.27)G′ (θ) は G(θ) を cos θ で微分したもので以下である.
[∂G(θ)
2c20 (1 + cos θ)G (θ) =
= a02
∂ cos θ[d20 + (1 + cos θ)2 ]′](2.28)

第2章

解析手法の基礎12

cos θ の位置ベクトルでの微分を考える. 原子 j, i, k の内角であるから
r ij · r ik
rij rik()()∂ cos θi1
cos θi r ij
1cos θi r ik
=
−
+
−
∂r i
rik
rij
rij
rij
rik
rik
cos θi =(2.29)(2.30)

次に Bji の微分を考える.

Bji = (1 + ζji )−δ より

∂Bji
∂ζji
= −δ (1 + ζji )−δ−1
∂r i
∂r i

(2.31)

中心∂ cos θj∂ζji j∑=fc (rjk )G′ (θj )
∂r i
∂r i
k̸=i,j(2.32)

r ji · r jk
rji rjk∂ cos θj
1 r jk cos θj r ji
=−
+
∂r i
rji rjkrji rjicos θj =(2.33)(2.34)

以上を整理すると[

F i = B¯ij VA′ (rij ) − VR′ (rij )

]r

ij

rij

(

)

θi
中心
(
) i∑
fc (rik )G′ (θi ) r1ik − cos
1
−δ−1
rij )
(
+ VA (rij ) −δ (1 + ζij )cos θi1′′2
k̸=i,j fc (rik )G(θi ) + fc (rik )G (θi ) rij − rik
r ji
θj
中心
(
) j∑
fc (rjk )G′ (θj ) cos
1
−δ−1
rji
rji
+ VA (rij ) −δ (1 + ζji )
r jk
1
′
2
−f
(r
)G
(θ
)
c jk
j rji
k̸=i,jrjk
原子 j の力はFj =[B¯ij VA′ (rij ) − VR′ (rij )

]

(

r ij
−
rij

)

[

1
∂Bij ∂Bji+ VA (rij )
+
2
∂r i
∂r i

r ij

]

(2.35)

である. Bij の微分の ζij について(中心
∂ cos θi∂ζij i∑=fc (rik )G′ (θi )
∂r i
∂r j
k̸=i,j)(2.36)∂ cos θi
1 r ik cos θi r ij
=−
+
∂r j
rij rikrij rij(中心
r jk∂ cos θj∂ζji j∑fc′ (rjk )G(θj )=+ fc (rjk )G′ (θj )
∂r j
rjk
∂r j
k̸=i,jrij

r ik(2.37))(2.38)rik第2章∂ cos θj
=
∂r j

(

1
cos θj−rjkrji

以上を整理すると

Fj =[B¯ij VA′ (rij ) − VR′ (rij )

]

(

r ij
−
rij

)

(

r ji
1
cos θj
+
−
rji
rji
rjk

解析手法の基礎
)

r jk
rjk

13

(2.39))

r ij

中心 f (r )G′ (θ ) cos θi
(
) i∑
1
c ik
i rij
−δ−1
+ VA (rij ) −δ (1 + ζij )1′2
k̸=i,j −fc (rik )G (θi ) rij

rij

r ik
rik

(

)

θj
r ji
fc (rjk )G′ (θj ) r1jk − cos
1
−δ−1
r
rji
ji
(
)
+ VA (rij ) −δ (1 + ζji )r jk
cos θj1′′2
k̸=i,j fc (rjk )G(θj ) + fc (rjk )G (θj ) rji − rjkrjk
原子 k の力は][1∂Bij ∂Bji+(2.40)
F k = VA (rij )
2
∂r k
∂r k

(

中心
) j∑

である. Bij , Bji の微分の ζij , ζji については
(

(

)

中心
∂ζij i∑r ik∂ cos θi=fc′ (rik )G(θi ) −+ fc (rik )G′ (θi )
∂r k
rik
∂r k
k̸=i,j

(

(

中心
r jk
∂ζji j∑fc′ (rjk )G(θj ) −
=
∂r k
rjk
k̸=i,j))∂ cos θj+ fc (rjk )G (θj )∂r k(2.41))′

である. 以上を整理すると
中心
(
) i∑
−fc (rik )G′((θi ) r1ik)
1
F k = VA (rij ) −δ (1 + ζij )−δ−1
cos θi′′2k̸=i,j −fc (rik )G(θi ) + fc (rik )G (θi )
rik
中心
(
) j∑
−fc (rjk )G′ (θj ) r1jk
1
−δ−1
+ VA (rij ) −δ (1 + ζji )cos θj′′2k̸=i,j −fc (rjk )G(θj ) + fc (rjk )G (θj ) rjk
同様に，Van der Waals ポテンシャルについては，

Fi = −∂Φ(rij )
∂r ij
∑



12
=−
4ϵ −
rij
j(̸=i)(

F j = −F i

σ
rij

)12

6
+
rij

(

σ
rij

)6 


r ij
rij

(2.42)

r ij
rij

r ik
rik
r ji
rji
r jk
rjk

(2.43)(2.44)

なお，添字 i，j ，k は原子をあらわし，rij , rik はそれぞれ粒子 i, j 間，i, k 間の距離，

θi は結合 rij , rik 間の角度を表している．

第2章

2.3

解析手法の基礎

14

高速化手法

原子数 N の系において粒子間の全相互作用を評価すると，1step に N × (N − 1) 回の計算が必要となり，N が大きくなると極めて膨大な計算量となる．実際には，一定距離以上離れた粒子は影響を及ぼさないので，作用を及ぼす範囲 (カットオフ半径 rc ) 内の粒子からの寄与を効率よく計算することにより高速化できる．従来よく用いられてきた高速化手法に粒子登録法がある．これは，図 2.4 に示したように，rc よりひとまわ
り大きい半径 rfc 内の粒子をメモリーに記憶し，その中で rc 内の相互作用を評価する方法であり，N × (rc 内粒子数 ≪ N − 1) に計算負荷が減少される．しかし，粒子登録法では rfc 半径より外の粒子が rc 内に達すると力の評価が適切でなくなるので，一定のステップ毎に登録粒子の更新 (N × (N − 1) 回の探査) を行わなければならない．このため，系がある程度の規模以上になると，粒子登録による高速化は登録更新の計算負荷により打ち消される．

rfc

rc

Fig.2.4 Schematic of bookkeeping method.

第2章

解析手法の基礎

15

別の高速化手法としてブロック分割法がある．図 2.5 に示すように，シミュレートする系をカットオフ距離程度の格子状に分割し，各ブロックに属する粒子をメモリーに記憶する．着目している粒子に作用する力を評価する際には，その粒子が属するブロックおよび隣接するブロックから相互作用する粒子を探索して行う．粒子が属する
ブロックは，粒子の位置座標をブロックの辺長 bx，by で除した際の整数により判断できるので，ブロック登録時の計算負荷は粒子数のオーダーとなる．したがって，粒子登録法では登録更新の負荷が大きくなるような大規模な系でも高速化が可能である．

y

by
0

x

bx

Fig.2.5 Schematic of domain decomposition method.

第2章

2.4

解析手法の基礎16

速度スケーリング法

分子動力学解析における温度制御には一般的には速度スケーリング法が用いられる．
この方法は，統計熱力学より導かれる式 (2.45) を用いて，以下のように制御する．

1 α α α 3m vi vi = k B T
22(2.45)

mα : 粒子αの質量
viα : 温度 T での粒子αの速度kB : Boltzmann 定数 = 1.38 × 10−23 [J/K](2.47)
目標の温度 T0 における原子 α の速度を viα0 とおくと viα0 は式 (2.48) のように表される．

(

viα0 =

3kB T0mα)0.5(2.48)

同様に，温度 T の時の原子 α の速度は式 (2.49) のように表される．
(

viα =

3kB T
mα

)0.5

(2.49)

よって，式 (2.48) と式 (2.49) より以下の式が得られる．
(

viα0
T0
=
α
vi
T

)0.5

(2.50)

つまり，系の温度を T から T0 にするには，式 (2.50) の右辺を現在の速度に掛けてやればよい．ただ，これだけでは原子配置に反映されないので，Verlet 法における ∆r αi (t+∆t)（式 2.51）を√

T0 /T ∆r αi (t + ∆t) と置き換える必要がある．

∆r αi (t + ∆t) = r αi (t + ∆t) − r αi (t) = r αi (t) − r αi (t − ∆t) + (∆t)2

F αi (t)(2.51)mα

平衡状態では，能勢の方法 [29] など外部との熱のやりとりをする変数を考慮した拡張系の分子動力学法によって得られるカノニカルアンサンブルに一致することが示されて
いる．

第2章

2.5

解析手法の基礎

17

カーボンナノチューブの幾何形状パラメータ

カーボンナノチューブの基本構造は六員環であるが，チューブ軸方向と六員環の並び方向の違いにより様々な形状が存在する．以後，チューブ形状を示すために，以下のような指標を用いる．周方向の六員環の個数が n1 ，一周するときに軸方向にずれる
六員環の個数が n2 であるチューブを (n1 ,n2 ) と表示する (図 2.6)．したがって，n2 ̸=0
はチューブが軸方向にらせん状になっていることを表す．

n2n1
Fig.2.6 Chiral index.

例えば図 2.7(a) の場合 (12, 0) ，図 2.7(b) の場合 (12, 4) と表される．

(a) (12,0)(b) (12,4)
Fig.2.7 Variety of carbon nanotubes.

第2章

解析手法の基礎18

n2 = 0 のときは zigzag 型，n1 > n2 (n2 ̸= 0) のときは chiral 型，n1 = n2 のときはarmchair 型と分類され，それぞれカイラリティにより性質が少し異なる特徴を持つ．
(n1 ,n2 ) チューブの半径 r は，原子間距離を a とすると以下の式で表される．
√r=

3 (n21 + n1 n2 + n22 ) a(2.52)

第3章
単層 CNC の構造安定性
3.1

小規模解析モデル

図 1.1 の SEM 像に示したように，CNC はコイル径が数十〜数百 [nm] のものが多く報告されている．しかし本研究では，多数のモデルで基礎的な検討を行うために，コ
イル径 100[nm] 以下の比較的小さな径で，巻き数が 1 の周期モデルについて解析を行う．短い CNT 切片を 36 角形状になるように角度をつけて接合箇所に欠陥を持たせて接合することで，CNC モデルを作成した．作成した CNC モデルの例を図 3.1 に示す．
カイラリティがらせん構造の安定性に与える影響を明確にするため，z 方向のみ周期境
界，xy 方向自由境界条件下で，カイラリティ(n1 ，n2 ) とチューブ半径 r の異なる 5 つのモデルについての検討を行う．作成したモデルの幾何形状パラメータを表 3.1 に示
す．コイル半径 R 及び周期長さ h は統一しており，Zigzag 型と Armchair 型について
は Chiral-2 と同程度のチューブ半径となるように作成した．

3.1.1

解析条件

初期構造緩和では，緩和前に熱処理を行うものと行わないものの二つの条件で計算を行った．熱処理を行わない場合では，温度を 10[K] として 100000[fs] の構造緩和を行った．熱処理を行う場合では，5000[fs] 毎に温度を 50[K] 上昇させ，500[K] に達した際に50000[fs] の間温度を保持し，その後上昇時と同様のペースで 10[K] となるまで温度を下降させた後，同様の 100000[fs] の構造緩和を行った．熱処理を行う際の温度制御を図3.2 に示す．z 方向に周期境界を与え，熱処理時には周期セルの長さを固定しているが，
19

第3章

単層 CNC の構造安定性

20100000[fs] の構造緩和時は z 方向の応力が 0 となるように z 方向セル長さを制御している．初期配置作成時の炭素原子間の結合距離は，今回使用した原子間ポテンシャルの平行粒子間距離である 0.145[nm] としている．各原子の初速度は Maxwell-Boltzmann
分布に従って乱数で与えている．温度制御は速度スケーリング法を用いており，制御の積分ステップは 1[fs] とした．

R
h

y

z

r

x

x
y

z

(a) top view(b) side viewFig.3.1 Example of SWCNC model.
Table 3.1 Parameter of small SWCNCs .
typechiral index

tube r [nm]coil R[nm]pitch h[nm]ZigzagArmchairChiral-1
Chiral-2Chiral-3(13，0)(8，8)(10，2)(10，5)(10，8)

0.52
0.55
0.45
0.53
0.62

5.0
5.0
5.0
5.0
5.0

3.0
3.0
3.0
3.0
3.0

Fig.3.2 Temperature control for heat treatment.

第3章

単層 CNC の構造安定性

21

解析結果

3.1.2

Zigzag-1 の構造緩和結果を図 3.3 に示す．(a) は熱処理なし，(b) は熱処理ありについての結果である．Zigzag-1 では，どちらにおいても安定したらせん構造を得ることは出来ず，z 方向応力 0 の下ではコイルが潰れてしまっていた．(a) 熱処理なしのチュー
ブの拡大図を見てみると，黒丸で示したチューブの接合箇所に四員環等の欠陥が存在してしまっていることが分かる．また，コイルの内側部分ではチューブに破断箇所が存在している．このように接合箇所での欠陥の影響が大きく，らせん構造が安定しなかったものと考えられる．(b) 熱処理ありの拡大図を見てみると，コイル内側での破断箇所は少し改善されているようであったが，黒丸の箇所に存在していた四員環などのチューブ構造に大きな変化は見られなかった．そのため，(b) 熱処理ありにおいてもらせん構造が安定しなかったものと考えられる．

z

y

x

x
y

z

(a)

z

y
x
z

x
y

(b)

Fig.3.3 Snapshots of Zigzag-1 after relaxation.

第3章

単層 CNC の構造安定性

22Armchair-1 の構造緩和結果を図 3.4 に示す．こちらも (a) 熱処理なし,(b) 熱処理あり
共に z 方向に潰れて，らせん構造を得ることは出来なかった．(a) 熱処理なしの拡大図を見てみると，Zigzag-1 と比べチューブに大きな破断は見られず，チューブ表面の凹
凸も少なかったが，四員環や大きな点欠陥は同様に存在していることが分かる．また，(b) 熱処理ありの拡大図と比較すると，四員環や点欠陥といった欠陥に大きな変化は見られなかった．

z

y
x

x
y

z

z

y
x
z

(a)

x
y

(b)
Fig.3.4 Snapshots of Armchair-1 after relaxation.

第3章

単層 CNC の構造安定性

23

Chiral-1 の構造緩和結果を図 3.5 に示す．このモデルも，共に安定したらせん構造を得ることは出来なかった．特に (b) 熱処理ありについては，チューブの半分程度まで破断している箇所がいくつか確認できた．また，二つの拡大図を比較しても構造に大きな変化は見られず，欠陥箇所の改善はされていなかった．このモデルは五つの内最もチューブ径の小さいモデルであるため，接合箇所の欠陥による影響を大きく受けてしまったものと考えられる．

z

y

x

x
y

z

(a)

z

y

x

x
z

y

(b)

Fig.3.5 Snapshots of Chiral-1 after relaxation.

第3章

単層 CNC の構造安定性

24

Chiral-2 の構造緩和結果を図 3.6 に示す．(a) 熱処理なしではらせん構造が潰れてしまっていたが，(b) 熱処理ありでは周期長さが 4.09[nm] と少し伸びた状態で潰れることなくらせん構造を得ることが出来た．また，Chiral-1 のようなチューブの破断が生じる
ことなく繋がっていることが分かる．(a) 熱処理なしの拡大図を見てみると点欠陥や四員環が残っておりチューブ表面の凹凸はあまり見られなかったのに対して，(b) 熱処理
ありの拡大図では，四員環は残っていたもの点欠陥が解消されており，またチューブ
表面に凹凸が生まれていた．また，(a) 熱処理なしではチューブが真っ直ぐ繋がっているのに対して，(b) 熱処理ありではチューブに矢印で示す方向にねじれが生じていた．

y

z

x

x
y

z

z

y
x
z

(a)

x
y

(b)
Fig.3.6 Snapshots of Chiral-2 after relaxation.

第3章

単層 CNC の構造安定性

25

Chiral-3 の構造緩和結果を図 3.7 に示す．このモデルでも (a) 熱処理なしでは安定したらせん構造が得られなかったのに対して，(b) 熱処理ありではほとんど周期長さが変わらない状態でらせん構造が安定していた．このモデルでも，(b) 熱処理ありの一部に大きな窪みは生じていたものチューブに破断は生じておらず，きれいに繋がっていることが確認できた．(a) 熱処理なしの拡大図を見てみると，このモデルにも点欠陥や四員環は存在していることが分かる．また，チューブに大きな凹凸は無く真っ直
ぐにチューブが繋がっていた．(b) 熱処理ありの拡大図を見ると，多少点欠陥が解消されており，またチューブ表面に凹凸が生じていた．また，Chiral-2(b) ほどでは無いがチューブに矢印の方向にねじれが生じている．

z

y
x

x

z

y

z

x
z

(a)

y

x
y

(b)
Fig.3.7 Snapshots of Chiral-3 after relaxation.

第3章

単層 CNC の構造安定性

26以上の結果から，まず，熱処理を与えない条件では安定したらせん構造を得ることが出来ないことが確認できた．また，500[K] での熱処理では点欠陥の改善はあったもの四員環の改善までは至らないことが分かった．Chiral-1 についてはチューブ径が小さく，接合箇所で破断が生じてしまっていためねじれについて判断することが出来な
いが，安定したらせん構造を得た Chiral-2(b) と Chiral-3(b) ではチューブにねじれが生じていたことが分かった．他の Zigzag，Armchair についてはねじれが生じておらず，真っ直ぐな状態で繋がっていたことを考えると安定したらせん構造を得るためには Chiral 型でありチューブにねじれを生じた状態でも安定出来ることが条件となるのではないかと考えられる．しかし，小規模のモデルではコイルに対して接合箇所の欠陥が与える影響が大きく，特に，Chiral-1 ではチューブが破断している箇所も存在していため，もう一回り規模の大きなモデルでの検討が必要と考える．

第3章

3.2

単層 CNC の構造安定性

27

大規模解析モデル

小規模のモデルでは，全体に対する接合箇所に含まれる欠陥構造の割合が大きく影
響が大きく出ていると考えられるため，規模を拡大したモデルを作成し，初期構造緩和シミュレーションを行った．作成したモデルの幾何形状パラメータを表 3.2 に示す．
先ほどの解析結果から小規模モデルでは Chiral 型でのみらせん構造が安定していたた
め，大規模モデルでは Chiral 型のみのモデルを作成した．3 種類の異なるカイラリティを持つモデルを作成し，周期長さを統一しコイル半径を変更した計 9 つのモデルについて解析を行う．

3.2.1

解析条件

先ほどの解析結果から熱処理を行った条件でのみ，らせん構造の安定したモデルが得られたことから，熱処理を与える構造緩和についてのみの解析とした．解析の条件
は先ほどと同様である．

Table 3.2 Parameter of Large SWCNCs .
typechiral index

Chiral-4Chiral-5Chiral-6Chiral-7
Chiral-8Chiral-9Chiral-10
Chiral-11
Chiral-12(30，20)(30，20)(30，20)(40，20)(40，20)(40，20)(50，30)(50，30)(50，30)tube r [nm] coil R[nm] pitch h[nm]
1.75
1.75
1.75
2.10
2.10
2.10
2.80
2.80
2.80

20.0
25.0
30.0
20.0
25.0
30.0
20.0
25.0
30.0

20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0

第3章

3.2.2

単層 CNC の構造安定性

28

解析結果

Chiral-4〜6 の構造緩和の結果を図 3.8 にまとめて示す．カイラリティが (30.20) のモ
デルでは，どのコイル径でもらせん構造が安定する結果となった．また，どのモデルでも緩和後に周期長さが大きくなっており，Chiral-4 は 25.8[nm]，Chiral-5 は 34.3[nm]，Chiral-6 は 42.0[nm] とコイル径が大きくなるにつれて周期長さも大きくなる傾向が見られた．図 3.9 に Chiral-4〜6 の一部を拡大して示す．どのモデルも接合箇所に点欠陥は見られたが，小規模モデルのように大きな凹凸やチューブに破断が生じることはな
かった．また，チューブの軸方向での断面はどのモデルでも楕円形になっており，前節の Chiral-2 のように接合箇所で大きくねじれているのではなく，チューブが少しずつ
全体でねじれていることが確認できた．接合箇所だけでねじれているのではなく，全体的に上方向にねじれが生じているため，コイル径の大きなモデルほどねじれによりz 方向周期長さが伸びたものだと考えられる．
Chiral-7〜9 の構造緩和の結果を図 3.10 に示す．Chiral-7 は，前節の小規模モデルのコイル径 5[nm] に比べると大きいが，らせん構造が安定せずに潰れている．Chiral-8,9については周期長さが 29.6[nm]，39.6[nm] とどちらも緩和前より伸びて，安定したらせん構造を保っている．また，コイル径が大きいほど z 方向周期長さ (ピッチ) が大きくなる傾向は先ほどの Chiral-4〜6 と同様である．図 3.11 に Chiral7〜9 の一部を拡大して示す．Chiral-7 では黒丸箇所でチューブに破断が生じており，チューブ自体にね
じれはほとんどなかった．それに対して Chiral-8,9 ではチューブに破断は起きておらず，Chiral-4〜6 と同様に少しずつ全体にねじれが起きていることが確認できた．
最後に，Chiral-10〜12 の構造緩和の結果を図 3.12 に示す．このカイラリティ(50，30)のものでも，コイル径が 20[nm] の Chiral-10 ではらせん構造が安定せず潰れてしまっていたが，Chiral11,12 ではそれぞれ周期長さが 25.4[nm]，30.8[nm] とやはり少し伸びた状態で安定していた．図 3.11 に Chiral7〜9 の一部を拡大して示す．Chiral-10 では，

Chiral-7 と同様に黒丸で示したコイルの一部の箇所に破断が生じており，チューブにねじれが生じていなかった．Chiral11,12 では接合箇所に破断は生じておらず，チュー
ブにねじれが起きる現象は同様であった．

第3章

単層 CNC の構造安定性

29以上の結果から，らせん構造が安定するモデルではチューブにねじれが生じている傾向が強いことが確認できた．また，チューブが上方向にねじれるためコイル径が大きいモデルほど，z 方向周期長さが大きい状態で安定する傾向があると推測される．小規模モデルの場合ではチューブの接合箇所で大きな凹凸が生じ，そこでねじれが起きていたのに対して，大規模モデルでは，接合箇所で大きな凹凸は見られず，ねじれも接合箇所に集中するのではなく，チューブが楕円形になるように全体でねじれていた．

第3章

単層 CNC の構造安定性

y
x
z

z
x
y

(a)Chiral-4(b)Chiral-5(c)Chiral-6

Fig.3.8 Snapshots of Chiral-4〜6 after relaxation.

(a)Chiral-4(b)Chiral-5(c)Chiral-6

Fig.3.9 Snapshots of zoom view of Chiral-4〜6 .

30

第3章

単層 CNC の構造安定性

y
x
z

z
x
y

(a)Chiral-7(b)Chiral-8(c)Chiral-9

Fig.3.10 Snapshots of Chiral-7〜9 after relaxation.

(a)Chiral-7(b)Chiral-8(c)Chiral-9

Fig.3.11 Snapshots of zoom view of Chiral-7〜9 .

31

第3章

単層 CNC の構造安定性

y
x
z

z
x
y

(a)Chiral-10(b)Chiral-11(c)Chiral-12

Fig.3.12 Snapshots of Chiral-10〜12 after relaxation.

(a)Chiral-10(b)Chiral-11(c)Chiral-12

Fig.3.13 Snapshots of zoom view of Chiral-10〜12 .

32

第4章
単層 CNC の引張シミュレーション
4.1

解析条件

前章の大規模モデルの中で，安定ならせん構造を得ることの出来なかった Chiral-7と 10 を除く 7 つのモデルを対象に，CNC の軸方向の周期長さを延伸することで引張シミュレーションを行う．7 つの解析モデルの構造緩和後の幾何形状パラメータを図4.1 に示す．引張はひずみ制御で行った．引張ひずみは毎ステップ増加させており，ひずみ速度を 1.0 × 10−5 [/fs] と 1.0 × 10−6 [/fs] の二つの条件で引張を行った．各原子の初速度は Maxwell Boltzmann 分布に従って乱数で与えており，温度は全て 10[K] で制御
し，積分の時間ステップは全て 1.0[fs] とした．

Table 4.1 Parameter of Large SWCNCs after relaxation.
typechiral index

Chiral-4Chiral-5Chiral-6Chiral-8Chiral-9Chiral-11
Chiral-12(30，20)(30，20)(30，20)(40，20)(40，20)(50，30)(50，30)tube r [nm] coil R[nm] pitch h[nm]
1.75
1.75
1.75
2.10
2.10
2.80
2.80

33

18.0
23.0
28.0
23.0
28.0
23.0
29.0

25.8
34.3
42.0
29.6
39.6
25.4
30.8

第4章

4.2

単層 CNC の引張シミュレーション

34

解析結果

図 4.1 に Chiral-4 の引張シミュレーション中の応力-ストレッチ関係を示す．(a) がひ
ずみ速度 1.0 × 10−5 [/fs] の結果，(b) がひずみ速度 1.0 × 10−6 [/fs] の結果である．また，
変形時の様子を図 4.2 と図 4.4 に示す．それぞれ，z 軸方向，y 軸方向及び，丸で示した箇所を拡大したものの 3 つである．どちらのひずみ速度でも λ=1.4 程度までは線形的に応力上昇しており，応力値はどちらもほとんど同じだった．また，この範囲でのばね定数は 8.94×10−4 [µN/nm] であった．その後，ひずみ速度の速い黒線では応力が大
きく揺らぎはじめ，λ=2.0 付近でピークを示し減少に転じる．図 4.2 を見ると，応力が揺らぎ始めた λ=1.5 付近では，赤丸で示した箇所で局所的な座屈を生じており，これらの発生が速い引張下での応力の揺らぎをもたらしたものと考えられる．また，図 4.3に図 4.2(b) に示したくぼみ箇所をさらに拡大したものを示す．点欠陥や五員環の箇所
では，チューブ表面に凹凸を生じていたが，局所的な座屈は五員環や点欠陥を含む接合箇所ではなく，赤線で示すように，六員環の並びに沿ってチューブ内側に折り込まれるような形で生じていることが確認できる．また，応力上昇が頭打ちとなる λ=2.0 では，黒丸で示した箇所にもねじれが集中し，著しく変形したため剛性が低下したものと考えられる．遅い引張では λ=1.4 付近で傾きが変わり，応力の揺らぎも比較的小さなまま λ=1.9 付近まで応力上昇しているが，やはりそこで頭打ちとなり減少に転じる．

λ=1.4〜1.9 でのばね定数は 4.69×10−4 [µN/nm] と 46%程度まで減少していた．応力上昇の鈍化した λ=1.4〜1.45 の間では，図 4.4(b) に赤丸で示したように局所的な座屈が発生している．また，応力が減少し始める λ=1.9 付近でも，特に黒丸の箇所で著しい変形集中が発生していた．しかし，どちらの速度でも λ=2.5 までの引張では，チュー
ブは破断しておらず，また，ひずみ速度の違いで応力の大まかな挙動や最大値，チュー
ブの窪んだ箇所に変化は見られなかった．
図 4.5〜図 4.8 に Chiral-5 の結果を示す．いずれも二次曲線的に応力上昇をしており，硬化していることがわかる．λ=1.4 までのばね定数は 3.39×10−4 [µN/nm]，λ=1.4 以降のばね定数は 1.19×10−3 [µN/nm] であり，およそ 3.5 倍にばね定数が上昇している．図第4章

単層 CNC の引張シミュレーション

35

4.7 に遅い引張での (a) 引張開始前，(b)λ=1.4，(c)λ=1.6 でのチューブの拡大図をそれぞれ示す．緩和前と比較すると，λ=1.4 では赤丸で示したチューブ接合箇所で z 軸方向に若干ではあるもののずれが生じている．しかし，λ=1.4 と λ=1.6 で比較するとチューブのずれに大きな違いは見られない．そのため，λ=1.4 までは引張によってチューブ
接合箇所でのずれが応力を解消したため，剛性は小さくなり，λ=1.4 以降はチューブのずれによる応力の解消は起こらず，チューブがチューブ軸方向に引張られるようになったため，剛性が上昇したものと考えられる．

λ=1.8 付近まではひずみ速度による応答の違いはほとんどない．図 4.6 と図 4.8 の
λ=1.8 における図を比較すると，遅い引張では λ=1.7〜1.8 で赤丸の箇所に局所的な座屈が生じていることが分かる．また，局所的な座屈が出来た箇所は Chiral-4 と同様に，チューブが内側に折り込まれるような形状をしており，結合が破断している箇所は見当たらない．
図 4.9〜図 4.11 に Chiral-6 の結果を示す．速い引張では大きな振動が現れているものの，λ=1.6〜1.7 近傍までは全体的な傾向は同じである．図 4.10(a) を見ると，λ=1.2と早い段階で局所的な座屈が赤丸の二か所にほぼ同時に発生しており，この影響により応力波が発生したものと考えられる．応力ピークは λ=1.8 近傍であるが，やはり
図 4.10(b) に黒丸で示した箇所で著しい変形が生じている．一方，遅い引張では，応力-ストレッチ曲線をよくみると初期は下に凸，後期は上に凸の曲線となっている．図4.11(a),(b) を見ると，下に凸すなわち硬化している λ=1.35〜1.4 で赤丸の箇所に局所的な座屈が生じている．λ=1.0〜1.65 でばね定数を算出するとひずみ速度に関わらず

1.09×10−3 [µN/nm] となる．
図 4.12〜図 4.14 に Chiral-8 のの結果を示す．これも λ=1.5 付近までは引張速度による差はない．λ=1.0〜1.45 でばね定数を算出すると 1.04×10−3 [µN/nm] である．速い引張では λ=1.5 で応力上昇が鈍化した後も下がらず，弾性的な応答を示している．一方，遅い引張では λ=1.9 以降応力が低下している．図 4.13 と図 4.14 で λ=1.5 以降の挙動を比較すると，速い引張では λ=2.1 でも黒丸のくびれた部分のチューブはまだ完
全に潰れていないが，遅い引張ではねじれが集中し，チューブが完全に潰れて折れている．見た目の構造は大きく変わらないので，構造の差というよりは変形速度による第4章

単層 CNC の引張シミュレーション

36

オーバーシュートと考えられる．
図 4.15〜図 4.17 に Chiral-9 の結果を示す．どちらも λ=1.4 付近までは速度による差がない．黒の変化は上に凸の二次曲線状で線形領域を決めるのが難しいが，λ=1.0 近傍の勾配からのばね定数は 1.20×10−3 [µN/nm] となる．速い引張では，ねじれが集中する箇所もあるが局所的な座屈が生じるのではなく，全体的にチューブがねじれによって潰れていく傾向が他のモデルより強く，楕円形に潰れていくことでねじれ剛性が小
さくなり上に凸の二次曲線的な挙動になったと考えられる．図 4.16 と図 4.17 を比較す
ると，y 軸方向からの図より，遅い引張では全体が一様に伸びているのではなく，一部にねじれが集中している箇所と直線的に伸びている箇所に分かれている．また，直線的に伸びている部分でのチューブ潰れの伝播が速かったため応力の挙動が急減したものと考えられる．
図 4.18〜図 4.20 に Chiral-11 の結果を示す．これも速度による差はあまりない．いずれも λ=1.9 付近でピークを生じ，その後直線的に減少しており「へ」の字の応答を示している．応力上昇時の勾配もこれまでのチューブに比べ小さく，ばね定数は8.32×10−4 [µN/nm] であった．このモデルが最も周期長さが小さいモデルであるため，引っ張られた距離が短かったこともばね定数が小さくなった要因であると考えられる．
図 4.19 と図 4.20 を見ると，どちらも他のモデルのように λ=1.4 から局所的な座屈が生じることはなく，チューブの潰れも小さかった．そのため，チューブのねじれ剛性の低下が抑えられ線形的に応力が上昇したと考えられる．
図 4.21〜図 4.23 に Chiral-12 の結果を示す．どちらも λ=1.9 付近までは速度による差はない．こちらも Chiral-11 と同様に応力上昇の勾配は小さく，λ=1.0〜1.9 でばね定数を算出すると 7.06×10−4 [µN/nm] である．速い引張では，λ=2.3 まで線形的な応力上昇を続け，大きな応力の減少はなかった．対して，遅い引張では λ=1.9 付近でピークを生じ，直線的に減少し，λ=2.1 でさらに応力が急減した．図 4.22 と図 4.23 を比較す
ると，遅い引張では λ=2.1 でチューブが潰れていたのに対して，速い引張ではチューブに座屈や潰れた箇所が生じていない．
以上のように，どのモデルについても引張初期の段階では，ひずみ速度による応力挙動への影響は小さい．二次曲線的な応力上昇を示すものもいくつかあるが，チューブ第4章

単層 CNC の引張シミュレーション

37

接合箇所がずれていくことによる応力緩和の影響やチューブが楕円状に変形すること
でねじれ剛性が小さくなる影響であると考えられる．応力が下がる原因としては，引張によるチューブのねじれによって，チューブに局所的な座屈が発生することやねじれによってチューブが潰れることが主な原因である．また，各モデルでの線形的な応力上昇領域でのばね定数を表 4.2 にまとめて示す．極端に値が大きいものや小さいものはなかったが，カイラリティ(50，30) のモデルが他と比較すると値が小さい傾向があった．しかし，これらのモデルは周期長さが統一されていないため，詳細な検討をするためには周期高さを統一したモデルの作成を行う必要があり，今後の課題である．

Table 4.2 Stiﬀness of SWCNCs .
typetube r [nm]coil R[nm]pitch h[nm]stiﬀness k [µN/nm]

Chiral-4Chiral-5Chiral-6Chiral-8Chiral-9Chiral-11
Chiral-12

1.75
1.75
1.75
2.10
2.10
2.80
2.80

18.0
23.0
28.0
23.0
28.0
23.0
29.0

25.8
34.3
42.0
29.6
39.6
25.4
30.8

8.94×10−4
1.19×10−3
1.09×10−3
1.04×10−3
1.20×10−3
8.32×10−4
7.06×10−4第4章

単層 CNC の引張シミュレーション

Stress , σ , GPa

0.15

0.1

0.05

(a)v =1.0×10-5(b)v =1.0×10-6

0

01

0.5
1.5

21

1.5
2.5

Stretch, λ , εRelativeelongation
Fig.4.1 Stress - stretch curves under tension (Chiral-4) .

38

第4章

単層 CNC の引張シミュレーション

z

y
x

x
y

z

(a)λ=1.5

y

z

x

x
y

z

(b)λ=1.6

z

y

x

x
z

y

(c)λ=2.0
Fig.4.2 Snapshots of Chiral-4 under tension (1.0 × 10−5 [/fs]) .

(a) Tube surface(b) Cross section

Fig.4.3 Zoom view of Chiral-4(b) (1.0 × 10−5 [/fs]) .

39

第4章

y

単層 CNC の引張シミュレーション

z

x
z

x
y

(a)λ=1.4

z

y
x

x
y

z

(b)λ=1.45

z

y
x

x
y

z

(c)λ=1.9

z

y
x
z

x
y

(d)λ=1.95
Fig.4.4 Snapshots of Chiral-4 under tension (1.0 × 10−6 [/fs]) .

40

第4章

単層 CNC の引張シミュレーション

Stress , σ , GPa

0.15

0.1

0.05

(a)v =1.0×10-5(b)v =1.0×10-6

0

01

0.5
1.5

12

1.5
2.5

Relativeelongation
Stretch, λ , εFig.4.5 Stress - stretch curves under tension (Chiral-5) .

41

第4章

y

単層 CNC の引張シミュレーション

z

x

x

z

y

(a)λ=1.8

y

z

x

x

z

y

(b)λ=2.0

y

z

x

x

z

y

(c)λ=2.1
Fig.4.6 Snapshots of Chiral-5 under tension (1.0 × 10−5 [/fs]) .

(a)λ=1.0(b)λ=1.4(c)λ=1.6

Fig.4.7 Zoom view of Chiral-5(b) (1.0 × 10−5 [/fs]) .

42

第4章

y

単層 CNC の引張シミュレーション

z

x

x

z

y

y

(a)λ=1.7

z

x

x

z

y

(b)λ=1.8

y

z

x
z

x
y

(c)λ=2.0
Fig.4.8 Snapshots of Chiral-5 under tension (1.0 × 10−6 [/fs]) .

43

第4章

単層 CNC の引張シミュレーション

Stress , σ , GPa

0.15

0.1

0.05

(a)v =1.0×10-5(b)v =1.0×10-6

0

01

0.5
1.5

21

1.5
2.5

Relativeelongation, εStretch, λFig.4.9 Stress - stretch curves under tension (Chiral-6) .

(i)

(ii)
(i)
y

z

x

x
y

z

(ii)

(a)λ=1.2

z

y
x
z

x
y

(b)λ=1.8
Fig.4.10 Snapshots of Chiral-6 under tension (1.0 × 10−5 [/fs]) .

44

第4章

y

単層 CNC の引張シミュレーション

z

x

x
y

z

(a)λ=1.35

z

y
x

x
y

z

(b)λ=1.4

z

y
x

x
y

z

(c)λ=1.6

y

z

x
z

x
y

(d)λ=1.9
Fig.4.11 Snapshots of Chiral-6 under tension (1.0 × 10−6 [/fs]) .

45

第4章

Stress , σ , GPa0.15

単層 CNC の引張シミュレーション(a)v =1.0×10-5(b)v =1.0×10-6

0.1

0.05

0

01

0.5
1.5

12

1.5
2.5

Stretch, λRelativeelongation, εFig.4.12 Stress - stretch curves under tension (Chiral-8) .

46

第4章

z

y
x

x
y

z

(a)λ=1.5

z

y
x

x
y

z

(b)λ=1.6zyxz

単層 CNC の引張シミュレーションxy(c)λ=2.1
Fig.4.13 Snapshots of Chiral-8 under tension (1.0 × 10−5 [/fs]) .

47

第4章

y

単層 CNC の引張シミュレーション

z

x

x
y

z

(a)λ=1.5

z

y
x

x
y

z

(b)λ=1.6

z

y
x

x
y

z

(c)λ=2.0

y

z

x
z

x
y

(d)λ=2.1
Fig.4.14 Snapshots of Chiral-8 under tension (1.0 × 10−6 [/fs]) .

48

第4章

単層 CNC の引張シミュレーション

Stress , σ , GPa

0.15

0.1

0.05

(a)v =1.0×10-5(b)v =1.0×10-6

0

01

0.5
1.5

12

1.5
2.5

Stretch, λRelativeelongation, εFig.4.15 Stress - stretch curves under tension (Chiral-9) .

49

第4章

単層 CNC の引張シミュレーション

z

y
x

x
y

z

y

(a)λ=1.4

z

x

x
y

z

(b)λ=1.5

y

z

x

x
y

z

(c)λ=2.1

z

y
x
z

x
y

(d)λ=2.2
Fig.4.16 Snapshots of Chiral-9 under tension (1.0 × 10−5 [/fs]) .

50

第4章

z

y
x

x
y

z

(a)λ=1.3

z

y
x

x
y

z

y

(b)λ=1.4zxz

単層 CNC の引張シミュレーションxy(c)λ=1.5
Fig.4.17 Snapshots of Chiral-9 under tension (1.0 × 10−6 [/fs]) .

51

第4章

単層 CNC の引張シミュレーション

Stress , σ , GPa

0.15(a)v =1.0×10-5(b)v =1.0×10-6
0.1

0.05

0

01

0.5
1.5

12

1.5
2.5

Relativeelongation, εStretch, λFig.4.18 Stress - stretch curves under tension (Chiral-11) .

52

第4章

z

y
x

x
y

z

(a)λ=1.4

z

y
x

x
y

z

y

(b)λ=1.5zxz

単層 CNC の引張シミュレーションxy(c)λ=2.0
Fig.4.19 Snapshots of Chiral-11 under tension (1.0 × 10−5 [/fs]) .

53

第4章

単層 CNC の引張シミュレーション

z

y
x

x
y

z

(a)λ=1.5

z

y
x

x
y

z

y

(b)λ=1.6

z

x

x
y

z

(c)λ=2.0

y

z

x
z

x
y

(d)λ=2.2
Fig.4.20 Snapshots of Chiral-11 under tension (1.0 × 10−6 [/fs]) .

54

第4章

Stress σ , GPa0.15

単層 CNC の引張シミュレーション(a)v =1.0×10-5(b)v =1.0×10-6

0.1

0.05

0

01

0.5
1.5

12

1.5
2.5

Stretch, λRelativeelongation, εFig.4.21 Stress - stretch curves under tension (Chiral-12) .

55

第4章

y

z

x

x
y

z

(a)λ=1.4

z

y
x

x
y

z

(b)λ=1.5zyxz

単層 CNC の引張シミュレーションxy(c)λ=2.1
Fig.4.22 Snapshots of Chiral-12 under tension (1.0 × 10−5 [/fs]) .

56

第4章

y

z

x

x
y

z

(a)λ=1.8

z

y
x

x
y

z

(b)λ=1.9zyxz

単層 CNC の引張シミュレーションxy(c)λ=2.1
Fig.4.23 Snapshots of Chiral-12 under tension (1.0 × 10−6 [/fs]) .

57

第5章
単層 CNC の繰り返し変形シミュ
レーション
5.1

解析条件

前章の引張の後に負荷ひずみ増分を反転させて元の周期長さまで戻す繰り返し変形
シミュレーションを行った．引張及び圧縮は前章と同様ひずみ制御で行っており，ひずみ速度を 1.0 × 10−5 [/fs] として計算を行った．各原子の初速度は Maxwell Boltzmann
分布に従って乱数で与えており，温度は全て 10[K] で制御し，積分の時間ステップは
全て 1.0[fs] とした．λ=2.0 まで引張を行い，その後同等のひずみ速度で λ=1.0 まで圧縮を行う条件とした．

5.2

解析結果

図 5.1 に Chiral-4 の繰り返し変形シミュレーション中の応力-ストレッチ関係を示す．
黒線が引張，赤線が除荷時の応力応答である．負荷反転時に大きく応力が揺らいでいた．チューブの構造に大きな変形等が生じている箇所は見当たらなかったことから，ひずみ速度が速いため，各原子にかかっている慣性の影響が大きく出たものだと考えられる．応力の揺らぎは λ=1.4 付近まで続いていたが，その後は，線形的な挙動になっており，黒線にほとんど重なっている．変形時の様子を図 5.2 に示す．赤丸の箇所に注目してみると，λ=2.0 の際にチューブに出来ていた局所的な座屈が，挙動が重なり始

58

第5章

単層 CNC の繰り返し変形シミュレーション

59

めた λ=1.4 まで圧縮しても残っていた．また，線形的に応力が減少している間に座屈は少しずつ解消されていたが，λ=1.0 でも座屈は完全には解消はされなかった．しかし，チューブの大部分は引張前と比較しても変化しておらず，チューブ変形の復元性が高いことが確認出来た．
図 5.3 と図 5.4 に Chiral-5 の結果を示す．Chiral-4 と比較すると負荷反転直後の応力の揺らぎは小さく，λ=1.7 付近から引張時の応力-ストレッチ曲線とほぼ同じになる．
図 5.4(b) と (c) を見ると，Chiral-4 と同様に，同ストレッチで引張にはなかった局所的な座屈が除荷時には存在していることが分かる．また，y 軸方向からの図より，除荷時の方がらせん形状に丸みを帯びているように見えることから，周期長さが長いため変形速度が速く，慣性から受ける影響は少し大きいと考えられる．しかし，Chiral-5 では λ=1.0 まで除荷すると先ほど存在していた座屈が完全に完全に解消されており，引張開始時とほぼ同じ形状に戻っている．先ほどの Chiral-4 でも局所的な座屈はおおよそ解消されていることからも，局所的な座屈に関しては復元性が高いことが再確認できる．
図 5.5 と図 5.6 に Chiral-6 の結果を示す．Chiral-6 では，負荷反転時に応力が大きく揺らいでおり，また，他モデルでは揺らぎながらも引張時の応力挙動に沿っていたのに対して，除荷時の応力が引張時よりも低い値となり，引張時応力-ストレッチの経路に戻ることはなかった．図 5.6 を見ると，λ=2.0 の際に赤丸の箇所で存在していた局
所的な座屈が除荷時に解消されていくのではなく，さらにねじれが大きくなっている．
このモデルは変形速度が速く，慣性の影響が強く働き，座屈の箇所が緩和されるよりも，座屈の箇所にさらにねじれが集中しまったためと考えられる．λ=1.0 まで除荷されると座屈箇所のねじれはある程度解消されているものの，全体的ならせん形状は他のモデルよりも乱れていた．
図 5.7 と図 5.8 に Chiral-8 の結果を示す．このモデルも負荷反転時は応力が揺らいでいたが，λ=1.4 付近から引張時の応力-ストレッチ曲線と重なった状態で線形的に応力が減少した．図 5.8 を見ると，λ=2.0 の時に発生していたねじれによる座屈が λ=1.4 付
近ですでに解消されており，チューブ表面に小さな凹凸が残っているのみであったた
め，他のモデルよりもより弾性的な応力挙動となった．

第5章

単層 CNC の繰り返し変形シミュレーション

60

図 5.9 と図 5.10 に Chiral-9 の結果を示す．やはり，周期長さが長いモデルでは，変形速度が速いため負荷反転時に慣性の影響を大きく受けている．引張時の応力-ストレッチ曲線と重なることはなく，また λ=1.0〜1.4 の除荷では応力値が負の値を示した．図5.10 を見ると，λ=2.0 の際に赤丸で示した箇所でチューブが潰れており，負荷反転時に，チューブの潰れが伝播したことがわかった．変形速度が速かったため，慣性の影
響も大きく，チューブの潰れの伝播はすぐにチューブ全体に伝わっていた．
図 5.11 と図 5.12 に Chiral-11 の結果を示す．周期長さが短く，変形速度も遅いため，
負荷反転時の慣性の影響は小さく，応力の揺れは小さかった．図 5.12 を見ると，λ=2.0の際もチューブにはねじれは生じているものの，局所的な座屈は起きておらず，除荷に切り替わるとすぐにねじれが解消されていた．元々弾性変形内での繰り返しである
ため，他のモデルよりもすぐに引張時の応力-ストレッチ曲線と重なっている．
図 5.13 と図 5.14 に Chiral-12 の結果を示す．こちらも負荷反転時に応力が振れただけで，ほとんど引張の応力-ストレッチ曲線と重なっている．弾性変形領域の場合，チュー
ブに生じるねじれも小さく，局所的な座屈もないため，負荷反転時の慣性の影響を除けば，引張と圧縮の応力-ストレッチ曲線はほとんど重なると考えられる．
以上のように，繰り返し変形シミュレーションの結果，引張時に生じる局所的な座屈は除荷時に解消され，チューブ変形の復元性が高いことがわかる．負荷反転時には応力が揺らいでいたが，負荷過程における弾性変形領域まで除荷すると応力-ストレッチ曲線は同じ経路に戻るものが多かった．しかし，引張の速度や量が大きく，チュー
ブ自体が潰れてしまう場合，圧縮時の応力挙動が引張の際より低い値を示すモデルも
あった．

第5章

単層 CNC の繰り返し変形シミュレーション

Stress , σ , GPa

0.2

0.1

0

Loading
Unloading

-0.1

-0.2

10

0.5
1.5

21

Stretch, λ , εRelativeelongation
Fig.5.1 Stress - stretch curves under tension and compression (Chiral-4) .

61

第5章

単層 CNC の繰り返し変形シミュレーション

z

y

x

x
y

z

(a)λ=1.0(Loading)

z

y

x

x
y

z

(b)λ=1.4(Loading)

z

y

x

x
y

z

(c)λ=2.0(Loading)

z

y

x

x
y

z

(d)λ=1.4(Unloading)

z

y

x

x
z

y

(e)λ=1.0(Unloading)Fig.5.2 Snapshots of Chiral-4 under tension and compression .

62

第5章

単層 CNC の繰り返し変形シミュレーション

Stress , σ , GPa

0.2

0.1

0

-0.1

-0.2

Loading
Unloading
10

0.5
1.5

21

Stretch, λRelativeelongation, εFig.5.3 Stress - stretch curves under tension and compression (Chiral-5) .

63

第5章

単層 CNC の繰り返し変形シミュレーション

z

y

x

x
y

z

(a)λ=1.0(Loading)

z

y

x

x
y

z

(b)λ=1.7(Loading)

z

y

x

x
y

z

(c)λ=1.7(Unloading)

z

y

x

x
z

y

(d)λ=1.0(Unloading)Fig.5.4 Snapshots of Chiral-5 under tension and compression .

64

第5章

単層 CNC の繰り返し変形シミュレーション

Stress , σ , GPa

0.2

0.1

0

Loading
Unloading

-0.1

-0.2

10

0.5
1.5

12

Relativeelongation, εStretch, λFig.5.5 Stress - stretch curves under tension and compression (Chiral-6) .

65

第5章

単層 CNC の繰り返し変形シミュレーション

z

y

x

x
y

z

(a)λ=1.0(Loading)

z

y

x

x
y

z

(b)λ=1.6(Loading)

z

y

x

x
y

z

(c)λ=2.0(Loading)

z

y

x

x
y

z

(d)λ=1.6(Unloading)

z

y

x

x
z

y

(e)λ=2.0(Unloading)Fig.5.6 Snapshots of Chiral-6 under tension and compression .

66

第5章

単層 CNC の繰り返し変形シミュレーション

Stress , σ , GPa

0.2

0.1

0

Loading
Unloading

-0.1

-0.2

10

0.5
1.5

12

Relatibe
elongation
Stretch, λ , εFig.5.7 Stress - stretch curves under tension and compression (Chiral-8) .

67

第5章

単層 CNC の繰り返し変形シミュレーション

z

y
x

x
y

z

(a)λ=1.0(Loading)

z

y

x

x
y

z

(b)λ=2.0(Loading)

z

y

x

x
y

z

(c)λ=1.4(Unloading)

z

y

x

x
z

y

(d)λ=1.0(Unloading)Fig.5.8 Snapshots of Chiral-8 under tension and compression .

68

第5章

単層 CNC の繰り返し変形シミュレーション

Stress , σ , GPa

0.2

0.1

0

Loading
Unloading

-0.1

-0.2

10

0.5
1.5

12

Relativeelongation, εStretch, λFig.5.9 Stress - stretch curves under tension and compression (Chiral-9) .

69

第5章

y

単層 CNC の繰り返し変形シミュレーション

z

x

x
y

z

(a)λ=1.0(Loading)

y

z

x

x
y

z

(b)λ=2.0(Loading)

z

y

x

x
y

z

(c)λ=1.6(Unloading)

z

y

x

x
z

y

(d)λ=1.0(Unloading)Fig.5.10 Snapshots of Chiral-9 under tension and compression .

70

第5章

単層 CNC の繰り返し変形シミュレーション

Stress , σ , GPa

0.2

0.1

0

Loading
Unloading

-0.1

-0.2

10

0.4
1.5

0.8
2

Stretch, λRelativeelongation, εFig.5.11 Stress - stretch curves under tension and compression (Chiral-11) .

71

第5章

y

単層 CNC の繰り返し変形シミュレーション

z

x

x
y

z

(a)λ=1.0(Loading)

y

z

x

x
y

z

(b)λ=2.0(Loading)

z

y
x
z

x
y

(d)λ=1.0(Unloading)Fig.5.12 Snapshots of Chiral-11 under tension and compression .

72

第5章

単層 CNC の繰り返し変形シミュレーション

Stress , σ , GPa

0.2

0.1

0

Loading
Unloading

-0.1

-0.2

10

0.5
1.5

21

Relativeelongation, εStretch, λFig.5.13 Stress - stretch curves under tension and compression (Chiral-12) .

73

第5章

単層 CNC の繰り返し変形シミュレーション

z

y

x

x
y

z

(a)λ=1.0(Loading)

z

y

x

x
y

z

(b)λ=2.0(Loading)

z

y
x
z

x
y

(c)λ=1.0(Unloading)Fig.5.14 Snapshots of Chiral-12 under tension and compression .

74

第6章

2 層 CNC の構造安定性・引張シミュ
レーション
6.1

初期構造緩和シミュレーション

6.1.1

解析条件

前章で用いた大規模モデルを重ねることで 2 層の MWCNC モデルを作成し，初期構造
緩和シミュレーションを行った．用いたモデルのパラメータを表 6.1 に示す．MWCNC-1
は前章の Chiral-5 と 8 を，MWCNC-2 は Chiral-5 と 11 を，MWCNC-3 は Chiral-8 と11 をそれぞれ重ねたモデルである．層間距離はそれぞれ 0.35[nm],1.15[nm],0.70[nm] となっており，MWCNC-1 についてはグラファイトの層間距離とほぼ同等とした．また，

MWCNC-4 は (40，20) との層間距離がグラファイトの層間距離と同等となるように(45，25) のチューブと重ねることで作成した．前章で，熱処理を与える条件でのみらせん構造が安定していたことから，今回も熱処理を与える条件での解析とした．熱処
理時の温度制御や緩和時間等は前章と同様の条件とした．

Table 6.1 Parameter of MWCNCs .
typechiral index

interlayer distance δ[nm]coil R[nm]pitch h[nm]

MWCNC-1
MWCNC-2
MWCNC-3
MWCNC-4(30，20), (40，20)(30，20), (50，30)(30，20), (50，30)(40，20), (45，25)

0.35
1.05
0.70
0.35

25.0
25.0
25.0
25.0

20.0
20.0
20.0
20.0

75

第6章

6.1.2

2 層 CNC の構造安定性・引張シミュレーション

76

解析結果

MWCNC-1 の構造緩和の結果を図 6.1 に示す．真ん中の拡大図 (a) は z 方向からの図の黒丸の箇所を拡大したもの，(b) はチューブ軸方向から見た断面図である．カイラ
リティ(30，20) と (40，20) を重ねたモデルでは，緩和後に z 方向周期長さが 33.9[nm]と少し伸びた状態でらせん構造が安定した．コイルの拡大図を見ると，MWCNC-1 はチューブ表面に大きな凹凸は生じておらず，また，ねじれが集中するような箇所は見られなかった．また，断面図を見ると，一定の層間距離を保った状態であり，赤丸で示した箇所に層間での結合が生じていることが確認出来る．結合が生まれた箇所は，各層での点欠陥や五員環が接近している箇所に多く見られ，ほとんどが接合部に存在していた．

MWCNC-2 の構造緩和の結果を図 6.2 に示す．カイラリティ(30，20) と (50，30) を重ねたモデルでは，チューブが全体的にねじれによって大きく変形しており，z 方向周
期長さが伸び続け，安定したらせん構造は得られなかった．特に，ねじれの集中していた黒丸の箇所では，赤丸に示すように内側と外側両方のチューブが完全に潰れている．チューブの層間距離が短い箇所で，層間での van der Waals の影響により，外側のチューブが内側のチューブを包み込むように変形していた．層間距離が 1.05[nm] と大きかったため，外側のチューブが変形が大きく，潰れてしまったことでらせん構造
も安定しなかったと考えられる．

MWCNC-3 の構造緩和の結果を図 6.3 に示す．カイラリティ(40，20) と (50，30) を重ねたモデルでも，チューブに潰れが生じ，構造緩和中 z 方向周期長さが伸び続けて
おり，安定したらせん構造を得ることは出来なかった．先ほどの MWCNC-2 と同様に黒丸の箇所にねじれが集中し，MWCNC-2 程ではないがチューブが潰れていることが確認出来る．MWCNC-3 は層間距離が 0.70[nm] と MWCNC-1 と比べるとやはり大きいため，MWCNC-2 と同様の原因で潰れたと考えられる．

MWCNC-4 の構造緩和の結果を図 6.4 に示す．カイラリティ(40，20) と (45，25) を重ねたモデルでは，チューブに潰れが生じることなく，z 方向周期長さが 28.0[nm] と少し伸びた状態で安定したらせん構造を得ることが出来た．MWCNC-1 と同様に，黒第6章

2 層 CNC の構造安定性・引張シミュレーション

77

丸の箇所を拡大してみてもチューブ表面に大きな凹凸やねじれの集中は見られなかった．一定の層間距離を保った状態で安定しており，赤丸の箇所でチューブの層間で結合が生じていることも確認できる．結合の生じた箇所は MWCNC-1 と同様に，点欠陥や五員環の接近している箇所であった．
以上の結果より，MWCNC-1 と 4 において安定したらせん構造を得ることが出来ていた．この二つのモデルは，層間距離をグラファイトの層間距離に近い 0.35[nm] としたモデルであり，層間距離が大きな MWCNC-2 と 3 においては，外側のチューブが内側のチューブを包み込むように潰れており，安定したらせん構造は得られなかったことから，層間距離の大きなものでは van der Waals の影響が大きく，安定したモデルを得ることが出来ないと考えられる．また，チューブ間で結合が生じている箇所が多数見られた．結合が生じた箇所は，チューブの接合箇所に存在する点欠陥や五員環の箇所で多く確認出来た．

第6章

2 層 CNC の構造安定性・引張シミュレーション(a) Tube surface

z

y
x

x

z

(b) Cross sectionyFig.6.1 Snapshots of MWCNC-1 after relaxation.

(a) Tube surface

y

z

x

x
y

z

(b) Cross section
Fig.6.2 Snapshots of MWCNC-2 after relaxation.

78

第6章

2 層 CNC の構造安定性・引張シミュレーション(a) Tube surface

y

z

x
z

x
y

(b) Cross section

Fig.6.3 Snapshots of MWCNC-3 after relaxation.

(a) Tube surface

y

z

x

x
y

z

(b) Cross section
Fig.6.4 Snapshots of MWCNC-4 after relaxation.

79

第6章

2 層 CNC の構造安定性・引張シミュレーション

6.2

引張シミュレーション

6.2.1

解析条件

80

前節で安定したらせん構造を得ることの出来た MWCNC-1 と 4 を対象に前章と同様の引張シミュレーションを行う．対象とした２つの解析モデルの構造緩和後のパラメー
タを表 6.2 に示す．ひずみ速度は 1.0 × 10−5 [/fs] のみとしており，その他の条件は第 4章の引張シミュレーションと同様とした．

Table 6.2 Parameter of MWCNCs after relaxation.
typechiral index

MWCNC-1 (30，20), (40，20)
MWCNC-4 (40，20), (45，25)

6.2.2

interlayer distance δ[nm]coil R[nm]pitch h[nm]

0.35
0.35

23.0
23.0

33.9
28.0

解析結果

図 6.5 に MWCNC-1 の引張シミュレーション中の応力-ストレッチ関係を示す．また，変形時の様子を図 6.6 に示す．λ=1.5 まで線形的に応力が上昇し，λ=1.5 付近で応力上昇が鈍化したものの上昇を続けた．λ=1.5 までのばね定数は 2.63×10−3 [µN/nm]，

λ=1.5 以降では 1.17×10−3 [µN/nm] であり，およそ 44%と半分以下まで減少した．カイラリティ(30，20) と (40，20) の単層モデルのばね定数はそれぞれ 1.19×10−3 [µN/nm]，

1.04×10−3 [µN/nm] であり，周期長さもそれほど違いがないが，初期応答は MWCNC-1の方がばね定数が大きい．これは層間での結合による相互作用によりねじれに対して強くなったためと考えられる．λ=1.8 程度で応力のピークを迎えた後，応力がわずかに減少し始め，λ=2.2 付近から応力が大きく減少した．図 6.7 に応力の鈍化した λ=1.5付近でのチューブ軸から見た断面図を示す．層間の結合に変化は見られず，応力が減少した λ=1.8 においても変化が見られなかった．応力のピーク値は，前章でのカイラ
リティ(30，20) の単層モデルと同程度の値であった．図 6.6 を見ると，応力上昇の勾配

第6章

2 層 CNC の構造安定性・引張シミュレーション

81が鈍化した λ=1.5 付近では，ねじれが集中し局所的な座屈が生じるような箇所はなく，黒丸の箇所を拡大しても，ねじれが少し強くなっているだけであった．λ=1.8 の時，赤丸で示した箇所で内側と外側両方のチューブに局所的な座屈が生じている．λ=2.2 では，黄丸で示した箇所にも，ねじれによる局所的な座屈が生じた．MWCNC でも，応力が減少する原因は SWCNC と同様にねじれの集中による局所的な座屈にあると考えられる．
図 6.8 と図 6.9 に MWCNC-4 の結果を示す．λ=1.6 までは線形的に応力が上昇しており，ばね定数を算出すると 2.65×10−3 [µN/nm] である．これは，MWCNC-1 とほぼ同じである．λ=1.6〜2.3 では応力挙動が横ばいとなっており，その後急減している．図6.9 を見ると，応力挙動が横ばいに変化した λ=1.5〜1.6 で赤丸の箇所でねじれが生じ始め，λ=2.0 ではそのねじれが大きくなっていることが分かる．λ=1.6〜2.3 での引張はこの箇所のねじれによって応力上昇しなくなり，横ばいの挙動になったものと考えられる．応力が減少に転じた λ=2.3 では，全体的にチューブが潰れ始めており，特に黒丸の箇所はチューブが完全に潰れ折れ曲がっていた．また，前章で扱った内側のカ
イラリティ(40.20) の単層チューブの応力挙動と応力挙動の形状は非常に似通っていたが応力のピークの値は MWCNC-4 の方が 2 倍程度大きな値となっていた．
以上の結果から，2 層の MWCNC においても局所的な座屈が生じることで応力が鈍化もしくは減少に転じることがわかる．また，ばね定数は複層化によって上昇し，その値は異なるカイラリティのチューブでほぼ等しい．ただし，今日はいずれのチュー
ブも (40，20) の同じカイラリティのチューブを有しているため，さらなる検討が必要
と考える．

第6章

2 層 CNC の構造安定性・引張シミュレーション

Stress , σ , GPa

0.15

0.1

0.05

0

10

0.5
1.5

21

1.5
2.5

23

Stretch, λRelativeelongation, εFig.6.5 Stress - stretch curves under tension (MWCNC-1) .

82

第6章

y

2 層 CNC の構造安定性・引張シミュレーション

z

x

x
y

z

(a)λ=1.4

y

z

x

x
y

z

(b)λ=1.5

y

z

x

x
y

z

(c)λ=1.8

y

z

x
z

x
y

(d)λ=2.2
Fig.6.6 Snapshots of MWCNC-1 under tension .

83

第6章

2 層 CNC の構造安定性・引張シミュレーション(a)λ=1.4(b)λ=1.5

Fig.6.7 Cross view of MWCNC-1 .

Stress , σ , GPa

0.15

0.1

0.05

0

10

0.5
1.5

21

1.5
2.5

23

Stretch, λRelativeelongation, εFig.6.8 Stress - stretch curves under tension (MWCNC-4) .

84

第6章

2 層 CNC の構造安定性・引張シミュレーション

z

y

x

x
y

z

(a)λ=1.5

z

y

x

x
y

z

(b)λ=1.6

z

y

x

x
y

z

(c)λ=2.0

z

y

x

x
z

y

(d)λ=2.3
Fig.6.9 Snapshots of MWCNC-4 under tension .

85

第7章
結論

本研究では，CNC の構造安定性ならびに力学特性について原子レベルからの知見を得るために，様々なカイラリティのストレートチューブ切片を接合して作成した単層の CNC について，分子動力学シミュレーションにより初期構造緩和，引張，繰り返し負荷などを行った．その後，2 層の複層 CNC についても，検討した．以下に，得られた結果を総括する．
第 2 章では，本研究で用いた解析手法の基礎について述べた．まず，分子動力学法の概要ならびに基礎方程式を示し，本研究で用いた数値積分法について説明した．次に，原子間相互作用の評価に用いられるポテンシャルエネルギーについて述べ，炭素原子に関するポテンシャル関数を具体的に説明した．さらに，大規模シミュレーショ
ンに必要な計算の高速化手法, ならびに CNT のカイラリティによる幾何形状について述べた．
第 3 章では，コイル径やカイラリティの異なる様々な単層 CNC を，コイル径 10[nm]の小規模な系と，コイル径 40〜60[nm] の大規模な径で構造緩和シミュレーションを行った．小規模な系では対象に 10[K] で構造緩和した場合と，500[K] まで昇温させてから 10[K] まで下げる熱処理をした場合について検討した．その結果，カイラリティのあるチューブに熱処理を与えた場合にのみ安定したらせん構造を得ることがわかった．Zigzag，Armchair で作成した CNC はコイル形状を保てなかった．安定した CNC
ではチューブがねじれながららせん形状を保つ傾向があった．熱処理を行わない条件
では，どのパターンにおいてもらせん構造の安定したモデルを得ることが出来なかっ

86

第7章

結論

87

た．この結果を受けて大規模モデルではカイラリティのあるチューブのみでかつ，熱処理条件下での検討とした．いくつかのチューブはコイル形状を保てなかったものの，多くは安定してコイル形状を保つことが出来た．周期長さ (コイルのピッチ) はコイル径が大きなものほど伸びており，チューブ全体に少しずつねじれを生じていた．
第 4 章では，第 3 章で安定した大規模 SWCNC を対象に引張シミュレーションを行った．ひずみ速度は 1.0 × 10−5 [/fs] と 1.0 × 10−6 [/fs] の二つの条件で行った．いずれの

CNC も，λ=1.5 程度までは線形的に応力上昇した．弾性変形領域では，振動には差があるものの応力上昇の勾配に違いはほとんど見られなかった．一部，応力が二次曲線的に上昇する硬化挙動を示す CNC も存在した．これは引張初期はチューブ接合部の変形により応力上昇が低くなったこと，引張後期はチューブ断面が潰れて楕円形状になり引張方向の剛性が変化したことなどが考えられる．他のチューブは応力上昇が鈍化し，ピーク応力を示して低下する．いずれのモデルも，最大応力を示した後の応力減少時には，チューブに生じるねじれによって座屈が生じていた．座屈は点欠陥や五員環を含む接合部ではなく，欠陥のない箇所で六員環の並びに沿って内側に折り込まれるよ
うな形で生じていた．引張初期のばね定数はいずれも 10−3 程度のオーダーであった．
第 5 章では，第 4 章の引張の後に負荷ひずみ増分を反転させて元の周期長さまで戻す繰り返し変形シミュレーションを行った．いずれも，負荷反転時に慣性により大きな応力振動を生じるが，負荷過程における弾性変形領域まで除荷すると応力-ストレッチ曲線は同じ経路に戻るものが多かった．これは引張時に生じていた局所的な座屈が除荷によって解消されるためで，元の周期長さまで除荷を終えると変形前の構造とほとんど違いがないものが多く，チューブ変形の復元性が高いことが確認出来た．しかし，周期長さが長い CNC では，変形速度が速く，負荷反転時の慣性による影響が大きい場合，局所的に座屈が生じている箇所にさらにねじれが集中し，チューブが完全に潰れてしまうケースも認められた．
第 6 章では，安定した大規模 SWCNC のうちから数種類の組み合わせで 2 層の MWCNCを作成し，構造緩和シミュレーション及び引張シミュレーションを行った．構造緩和の結果，チューブ間の層間距離がグラファイトの層間隔 0.35[nm] より大きい (外のチュー
ブとうちのチューブの半径差が大きい)MWCNC では，内側のチューブを外側のチュー

第7章

結論

88

ブが包み込むように接近し，安定なコイル形状を保てなかった．一方，層間距離がグラファイトの層間隔に近い MWCNC では，構造緩和後も一定の層間距離を保ったまま安定したらせん構造を得ることが出来ていた．ただし，チューブの層間では一部結合を生じており，点欠陥や五員環を含む接合箇所に多く見られた．コイル径 46[nm]，カイラ
リティ(30，20)-(40，20) と (40，20)-(40，25) の MWCNC を引っ張った結果，いずれもばね定数は 2.65×10−3 [µN/nm] 程度となり，SWCNC のばね定数よりも大きくなった．また，MWCNC においてもピーク応力後の応力減少はチューブにねじれによって局所的な座屈を生じることによってもたらされていた．

参考文献[1] Kroto, H.W., Heath,J.R., O’Brien, S.C., Curl, R.F., and Smalley, R.E., Nature ,318, (1985), 162.
[2] Iijima, S., Nature, 354, (1991), 56.
[3] 中山, 表面科学, 6, (2004), 332.
[4] 須田, 田上, 滝川, J. Plasma Fusion Res, 88, (2012), 629.
[5] Pan, L., Zhang, M., and Nakayama, Y., J, Appl.Phts., 91, (2002), 10058.
[6] Akagi, K., Tamura, R., and Tsukada, M., Phys.Rev.Lett., 74, (1995), 2307.
[7] Davis, W.R., Slawson, R.J., and Rigby, G.R., Nature, 171, (1953), 756.
[8] R.T.K. Baker. et al, Nature, 253, (1975), 379.
[9] S. Motojima et al, Appl. Phys. Lett, 56, (1990), 321.
[10] S. Iijima, and T. Ichihashi., Nature, 354, (1993), 603.
[11] A.G. Nasibulin et al, Carbon, 43, (2005), 2251.
[12] Kawakubo, Y., Gotou, M., Toumatsu, M., and Miyazaki, K., Proceedings of
WTC2005, Washington, USA, (2005), WTC2005-53602.
[13] T. Kodama et al, Abstract of the 41st Fullerenes Nanotubes Graphene General
Symposium, (2011), 7.
[14] 小山, 繊維学会誌, 62, (2006), 372.
[15] 尾方, 花生, 渋谷, 材料, 55, (2006), 754.
[16] 渋田, 丸山, 日本機械学会論文集 (B 編), 72, (2006), 2513.

89

参考文献

90

[17] Hirai, Y. et al, Jpn J Appl Phys, 42 6B, (2003), 4120.
[18] Deguchi, H., Yamaguchi, Y. et al, Chem. Phys. Lett., 503, (2011), 272-276.
[19] 西村, チュン, 荒井, 計算数理工学論文集, 10, (2010), 13.
[20] Pan, L., Zhang, M., and Nakayama, Y., J, Appl. Phys., 91, (2002), 10058.
[21] 細川, 後藤, 李, 長坂, 土屋, 末金, 野坂, 潘, 中山, 第 50 回応用物理学関係連合関連
会講演予稿集, 3, (2003), 601.

[22] Nakayama, Y., Pan, L., Suekane, O., and Nosaka, T., Proc. Of IS and T’s NIP
19 Int. Conf. on Digital Printing Tech., (2003), 727.
[23] T. Katsumata et al., Trans. Mater. Res. Soc. Jpn., 29, (2004), 501.
[24] 久米, 長谷川, 野坂, 中山, J. Jpn. Soc. Powder Metallurgy, 58, (2010), 3.
[25] 香川, 山口, 平原, 計算力学講演会講演論文集, 23, (2010), 575.
[26] Lizhao, Liu., Jijun, Zhao., Syntheses and Applications of Carbon Nanotubes andTheir Composites, (2013), 257.
[27] Brenner, D.W., Phys.Rev.B., 42, 15, (1990), 9458.
[28] Ahmadieh, A.A., Raﬁzadeh, H.A., Phys.Rev.B., 7, (1973), 4527.
[29] 上田顯, コンピューターシミュレーション, (1990), 朝倉書店.

学術論文・学術講演
学術講演

 分子動力学法による MW カーボンナノコイルの作成と引張シミュレーション
毛利

友宙，西村

英晃，屋代

如月

日本材料学会 第 18 回分子動力学シンポジウム, 東京工業大学, (2013.5)

91学術論文・学術講演
論文貼り付け用ページ

92学術論文・学術講演
論文貼り付け用ページ93

謝辞
本研究を遂行するにあたり，仔細に渡る懇切丁寧な御指導を賜りました屋代如月教
授 (現 岐阜大学) に，心より感謝の意を表します．浅学非才の著者に対し，三年間に渡り親切に，時に厳しく御指導頂いたことは人生の大きな財産となりました．また，本論文を完成するにあたり，研究報告会などにおいて様々な御助言を賜りました田中克
志教授，日々の研究活動が円滑に行えるよう便宜を図っていただきました片山雷太技
術職員に深く感謝いたします．
若輩だった私に貴重な御助言を頂き，有益な議論を交わした西村英晃氏 (現 阪神高速道路株式会社) に心より御礼申し上げます．日々の学生生活や大学院一年次の就職活動において，多くの有益なアドバイスを与えてくださった表面・界面工学研究室の諸先輩方に厚く御礼申し上げます．
同じ研究室に配属されて以来共に研究に取り組み，互いに切磋琢磨し，時には励ましあった伊藤貴彦君，小出卓矢君，津川悠太君，菱谷卓也君，藤麻成貴君，本馬剛徳
君，森太一君に御礼申し上げます．また，同じ材料系として議論を交わし，多くの時間を共有した固体力学研究室の同輩諸子に御礼申し上げます．
最後に，六年間の学生生活を暖かく見守り，経済的にも精神的にも支えて頂いた家族に最大の敬意と感謝の意を表します．ありがとうございました.
平成 27 年 2 月
毛利94

友宙