修士論文

シリコンにおける非弾性変形開始時の局所力学状態に関する分子動力学的研究

指導教員：屋代 如月

大石 直樹

2010 年 2 月
神戸大学大学院 工学研究科 博士課程前期課程 機械工学専攻

Molecular Dynamics Study on
Local Mechanical Condition at the Onset ofNon-elastic Deformation in Silicon

February 2010
Department of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, Japan
Naoki OISHI

要約
本論文では，シリコンの変形・破壊メカニズムについて原子レベルから知見を得るため，様々な形状のシリコンについて分子動力学法による変形シミュレーションを行い，非弾性変形挙動発生時の局所安定性について議論した．まず，用いた Tersoﬀ 型ポテン
シャルについて，理想完全結晶状態で [001] 方向に単軸引張を与えたときの格子不安定条件を静力学解析により求めた．横ひずみ ε´=0，横応力 σ´=0 のいずれにおいても最大応力を示すひずみより前に格子不安定となるひずみが存在した．次に分子動力学シ
ミュレーションによりバルクおよびナノワイヤでの [001] 方向への引張を行った．バルクは内部にループ状の欠陥を生成して応力急減したが，各原子の Bijα の正値性を調べた結果，欠陥生成直前に不安定原子が発生していた．ナノワイヤの引張では、4 つあるエッジのうち，後から不安定原子が発生したエッジ部に転位が発生していた．その理由として，エッジ部を「引き止める」Si 原子列があるか否かによることを明らかにした．次に，無限平板およびナノワイヤの 3 点曲げシミュレーションを行った．無限平板，ナノワイヤいずれも引張を受ける表面からのぜい性的な破壊を生じたが，小さい
無限平板では引張側表面の原子構造が乱れただけだった．全原子の Bijα を調べたところ，最終破断部よりむしろ単純支持点近傍で非常に多くの原子が不安定となっていた．
そこで，主応力成分に相当する 3×3 行列の正値性のみを見たところ，最終破断を生じる領域に不安定原子群が認められた．最後に軸方向の結晶方位が異なる単結晶ならびに bamboo 状のカンチレバーの曲げシミュレーションを行ったが，対象とした寸法な
らびに曲げ条件では非弾性変形を生じなかった．3×3 主小行列式の値で不安定評価をしても当然ながら大きな不安定領域は現れていない．また，結晶方位を [100][010][001]から回転させた座標系では，Bijα の固有値解析が必要であることも示した．

Summary
For a new insight on deformation and fracture mechanism of silicon, various shaped
Si are deformed in the molecular dynamics simulations, and the onset condition of
inelastic deformation is discussed based on the ``local lattice stability´´. First, theelastic limit of Si under the [001] uniaxial tension is calculated by statics analysis or
the convensional lattice stability analysis for Tersoﬀ potential. The instability point
emerges before the stress-strain peak for both the transverse conditions of ε´=0 andσ´=0. Then bulk and nano wire of silicon are subjected to the [001] unixial tention
by molecular dynamics simulation. Loop-like defects nucleate in the bulk silicon andlead the stress drop (inelastic deformation). We have evaluated the atomic stability
from the positiveness of the atomic elastic stiﬀness and revealed that unstable atoms
emerges just before the stress drop. In the simulation of Si nano-wire, dislocations
generated from two edges of four surfaces. It is noteworthy that the unstable atoms
are found on the other two edges faster and more than the deformed edges. Detailobservation reveals that the edges have atomic line just on the edge while the deformed
ones have rounded shape. Next, three-point-bending simulations are performed on an
infinite plate and nano bar. Both the plate and bar show sudden cracking by brittle
deformation from the tensile surface under the center punch. However, if we use asmall simulation cell for infinite plate the unstable behavior is stopped only with thedisorder of the tensile surface. When Bijα for all the atoms is examined, a lot of atoms
are judged as ``unstable´´ around the simple support of 3-point bending, although theunstable fracture starts beneath the center punch. Thus we evaluated the 3×3 minor
determinant in the 6×6 Bijα matrix, corresponding to the principal stress direction.
We can find unstable atoms at final breaking point with this ciriteria. Finaly, thebending simulation of nano cantilevers are implemented. The cantilever is made fromsingle crystal Si or multi grain bamboo structure. For the single crystal cantilever, thecrystal orientation is rotated to [100], [110] and [111] for the cantilever axis. However,
all the cantilever does not show inelastic deformation and there is little diﬀerence inthe force-displacement responses. We have checked the positiveness of 3×3 matrix ofBijα , however, needless to say there is no domain of unstable atoms since we couldn’t
observe the inelastic deformation. Here, it is also pointed out that we should use
the eigenvalue analysis for Bijα , if the crystallographic orientation is changed from theprincipal stress axis.

目次
第1章 緒 論1

第 2 章 解析手法

42.1

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

4

2.2

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

5

2.3

Tersoﬀ 型ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.4

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

10

2.5

高速化手法

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.6

格子不安定性解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.6.1

Born の安定基準 . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.6.2

ブラベー格子を用いた安定性解析 . . . . . . . . . . . . . . . . .

13

2.6.3

弾性剛性係数による安定性解析 . . . . . . . . . . . . . . . . . .

13

第 3 章 [001] 引張シミュレーションならびに格子不安定解析

3.1

3.2

静力学解析

16
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

3.1.1

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

16

3.1.2

応力−ひずみ関係 . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.1.3

格子安定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

分子動力学シミュレーション. . . . . . . . . . . . . . . . . . . . . . .

21

3.2.1

シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . .

21

3.2.2

シミュレーション結果と考察 . . . . . . . . . . . . . . . . . . .

22

3.2.3

局所格子不安定性 . . . . . . . . . . . . . . . . . . . . . . . . . .

27

第 4 章 曲げシミュレーションならびにi目次
格子不安定解析

4.1

4.2

ii
31

シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

4.1.1

無限平板モデル . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

4.1.2

ナノワイヤモデル . . . . . . . . . . . . . . . . . . . . . . . . . .

33

シミュレーション結果と考察. . . . . . . . . . . . . . . . . . . . . . .

34

4.2.1

無限平板の変形挙動 . . . . . . . . . . . . . . . . . . . . . . . .

34

4.2.2

ナノワイヤの変形挙動 . . . . . . . . . . . . . . . . . . . . . . .

43

4.2.3

局所格子不安定性 . . . . . . . . . . . . . . . . . . . . . . . . . .

45

第 5 章 ナノカンチレバーの曲げシミュレーションならびに格子不安定解析

53

5.1

シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

5.2

シミュレーション結果と考察. . . . . . . . . . . . . . . . . . . . . . .

55

5.2.1

押し込み荷重−変位関係 . . . . . . . . . . . . . . . . . . . . . .

55

5.2.2

局所格子不安定性 . . . . . . . . . . . . . . . . . . . . . . . . . .

62

第 6 章 結論

65

参考文献

68

関連発表論文・講演論文

73

謝辞

77

第1章
緒 論

半導体産業の発展は，微細加工技術の進歩と切り離すことが出来ない．1980 年代後半に，半導体微細加工技術を利用したモーターや歯車列などの微細な機械的部品の製作技術が米国で発表されて以来，世界各国でマイクロマシン技術に関する研究が活発に行われている．特に 3 次元的微細化の代表的なものとして機械要素部品，センサー，アクチュエータ，電子回路を一つの基板上に集積化した MEMS(Micro Electro Mechanical

Systems) が挙げられる．この MEMS 技術を利用してマイクロモーターや加速度セン
サー，カプセル内視鏡など広範囲にわたるデバイスの作製が可能となっている．
単結晶シリコンは MEMS において最も共通な構成物質の一つであり，マイクロ・ナ
ノスケールでの材料特性の研究が活発に行われている．原子間力顕微鏡を利用した単結晶シリコンナノワイヤの曲げ試験によって，低温では弾性変形から急にぜい性的な破壊が生じ，高温になるにつれ破断前に塑性変形する領域が大きくなること，常温に近い温度で転位が発生することが報告されている (1) ．シリコンマイクロカンチレバー
押し込み試験では，カンチレバーの寸法は応力−ひずみ関係だけでなく，曲げ強度やヤング率等にも影響すること (2) が報告されている．他にも，単結晶シリコン薄膜の単軸引張試験による結晶方位の破壊への影響 (3) ，多結晶シリコン薄膜の引張強度評価 (4)などが行われている．
また，近年の計算機能力の向上を背景に，原子シミュレーションによる変形・破壊現象の解明も盛んに行われている．Jing らは分子動力学法を用いた単結晶シリコンナノ
ワイヤ (SiNWs) の単軸引張，単軸圧縮シミュレーションを行っている (5), (6) ．単軸引張12シミュレーションでは臨界荷重は温度の増加またひずみ率の増加で小さくなり，ナノワイヤの直径の増加によって大きくなること，軸方向の結晶方位が異なる単結晶ナノワイヤの単軸圧縮シミュレーションでは結晶方位によって変形挙動が異なること，などを報告している．Park らはアモルファスシリコンの様々な温度で結晶化過程をシミュ
レーション (8) するとともに，シリコンナノカンチレバーの曲げと水平振動によって弾性特性のサイズ効果を検討している (7) ．Izumi らはシリコン基板の原子ステップからの 3 次元転位ループの核生成について，活性化エネルギーとそのエネルギーバリアのsaddle-point 形状を評価している (9) ．また，藤井らはき裂を有する単結晶シリコン薄膜の単軸引張シミュレーションを行い，切欠き深さが増加するほど弾性率および最大引張応力が低下すること (11) を報告している．なお，対象とした原子数が少ないため
直接変形・破壊を扱ったものではないが，Si の基礎的物性の評価として第一原理計算による検討も広くなされており，単結晶シリコンの理想せん断強度 (12), (13) や理想引張強度等 (14) が多数報告されている．
本研究では，バルク，ナノワイヤ，無限平板そしてカンチレバーなど，様々な形状のシリコンについて単軸引張や曲げ変形シミュレーションを分子動力学法を用いて行うとともに，非弾性変形開始時の局所格子不安定性 (15) について検討する．局所格子安定性とは，各原子位置におけるエネルギーの空間勾配に相当する物理量 (原子弾性剛性係数 Bijα ) の正値性を指し，これまでの研究により detBijα ＜ 0 の正負で単結晶表面からの転位発生 (15), (17) や微視的へき開破壊の発生 (16) ，ナノインデンテーションにおける圧子下の転位発生 (18) などが評価できることが示されている．以下，各章の概略を示す．
第 2 章では解析手法について述べる．はじめに，分子動力学法の基礎式を示し，本研究で原子間相互作用の評価に用いた Tersoﬀ 型ポテンシャルの概要について説明する．また，大規模分子動力学計算の手法である領域分割による高速化について説明する．さらに，各原子の安定性を評価するために用いた局所格子不安定性解析について概説する．第 3 章では，Tersoﬀ 型ポテンシャル (21), (22) を用いたときの単結晶シリコンの [001] 方向引張における格子不安定条件を静力学解析で求めるとともに，分子動力学法で [001] 方向に単軸引張りシミュレーションを行い，非弾性変形発生時の局所格子不3
安定性について検討する．第 4 章では，単結晶シリコンの無限平板およびナノワイヤの曲げシミュレーションを行い，寸法または温度による変形挙動の違いと，非弾性変形発生時の局所格子不安定性について検討する．第 5 章では，単結晶で異なる結晶方位としたカンチレバーならびに bamboo 状に結晶を有するカンチレバーの押し込みシ
ミュレーションを行い，その変形挙動と局所格子不安定性について検討する．最後に，第 6 章で本研究の総括を述べる．

第2章
解析手法

本章では，分子動力学法の基礎方程式を示し，本研究で原子間相互作用効果の評価に用いた Tersoﬀ 型ポテンシャルの概要について説明する．続いて分子動力学計算における温度制御のために用いた速度スケーリング法，領域分割による高速化手法について述べる．最後に，局所の安定性を評価するために用いた局所格子不安定解析について概説する．

2.1

分子動力学法

分子動力学法 (Molecular Dynamics Method ; MD) は，系を構成する個々の原子についてニュートンの運動方程式mαd2 r α
= Fα
2
dt

(2.1)

を作成し，これを数値積分することによって全原子の運動を追跡する手法である．ここで t は時間，r α ，mα はそれぞれ原子 α の位置ベクトル及び質量である．原子 α に作用する力 F α は系全体のポテンシャルエネルギー Etot の空間座標についての勾配ベクトルから次式のように求められる．

Fα = −

∂Etot
∂r α

(2.2)

式 (2.1) の数値積分には，Verlet の方法が簡便で高精度が得られるため MD 法ではよく用いられる．時刻 t±∆t での原子 α の座標 r α (t ± ∆t) を Taylor 展開すると

r α (t ± ∆t) = r α (t) ± ∆t

dr α (∆t)2 d2 r α (∆t)3 d3 r α+±+ O((∆t)4 )dt
2 dt2
3! dt34(2.3)

2.2 原子間ポテンシャル

5となる．両式の和をとり式 (2.1) を代入すると
2F

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

α

α

α

(t)
+ O((∆t)4 )
m

(2.4)

を得る．これより，時刻 t − ∆t と t における全原子の位置を既知として，時刻 t + ∆t
における任意の原子 α の位置を求めることができる．

2.2

原子間ポテンシャル

系のポテンシャルエネルギー Etot の評価は，以下の 3 つに大別される．

(1) 経験的ポテンシャル(2) 半経験的ポテンシャル(3) 非経験的手法 (第一原理計算)
経験的ポテンシャルは，量子力学の厳密な理論に基づいて決定されるのではなく，ポテンシャルを微分可能な未定係数を含む簡単な関数形で仮定し，従来の実験的事実に合致するようにその未定係数が決められる．半経験的ポテンシャルは，密度汎関数論
より導出される形で定義されるが，そのポテンシャルパラメータは平衡状態でのマク
ロな物性値や，あるいは ab-initio な計算により求められた値に対してフィッティング
される．非経験的手法とは，従来の特性値などを一切用いず，原子核の位置ならびに種類のみを必要情報とし，各時刻における電子状態を量子力学に基づいて解くことで，遂次原子に働く力を精密に評価する手法である．

MD 法において，原子 α に作用する力 F α は系のエネルギー Etot の空間微分によって求めるため (式 (2.2))，系のポテンシャルエネルギー Etot をいかに精度よく評価するかが重要となる．(3) の第一原理分子動力学法は，計算量が極めて膨大になるため，変形・破壊のような多数の原子の動的挙動への直接的な適用は困難である．そこで，原子間相互作用を簡略評価する (1)，(2) の原子間ポテンシャルが通常用いられる．

2.3 Tersoﬀ 型ポテンシャル

2.3

6

Tersoﬀ 型ポテンシャル

Si や C 等の共有結合材料は，sp3 の強い方向性を持った結合を有するため単純な 2 体
間ポテンシャルでは表現できない．そこで結合角の効果を取り入れた多体ポテンシャ
ルを用いる必要がある．共有結合材料の代表的なポテンシャルである Tersoﬀ 型ポテン
シャル (21)–(23) は，経験的ボンドオーダーポテンシャルと呼ばれ，配位数や結合角などの周りの環境に合わせて結合状態を変化させる．また，Si のポテンシャルパラメータには，多形態のシリコンの凝集エネルギー等がフィッティングされており，様々な条件下での原子構造を柔軟に表現できるようにしている．

Tersoﬀ 型ポテンシャルでは，系全体の結合エネルギーは各原子ごとのエネルギーの和として次のように表わされる．
N ∑
N
1∑
E αβ
Etot =E =2 α β(̸=α)
α
N
∑

α

(2.5)

各原子ごとのエネルギー E αβ は次のように表わされる．
[

(2.6)

fR = A exp(−λ1 rαβ )(2.7)

fA = −B exp(−λ2 rαβ )

(2.8)

1
(bαβ )∗ = (bαβ + bβα )
2

(2.9)



1,


) 

 (


αβ

π
r
−
R
 /2
fc (rαβ ) = 1 − sin 

2D






0,

]

E αβ = fc (rαβ ) fR (rαβ ) + (bαβ )∗ · fA (rαβ ): rαβ < R − D: R − D < rαβ < R + D(2.10): rαβ > R + D

ここで fR (rαβ ) は斥力を表わす項であり，(bαβ )∗ · fA (rαβ ) は引力を表わす項である．
また fc (rαβ ) はカットオフを表わす項である．bαβ は，原子 α，原子 β 以外に原子 γ も含めた 3 個の原子によって定められる項であり，次式で表わされる．

bαβ = {1 + (β α )n (ζ αβ )n }−1/2n(2.11)

2.3 Tersoﬀ 型ポテンシャル
ζ αβ =∑[

γ(̸=α,β)

g(θ

αβγ

(

fc (rαγ )g(θαβγ )exp λ33 rαβ − rαγ[() = 1 + c /d − c / d + h − cos θ
2

2

2

2

)]

)
αβγ 2

7

(2.12)](2.13)

ここで θαβγ は結合 rαβ ，rαγ 間の角度を示している (図 2.1)．

γ

θαβγ
α

β

Fig.2.1 Three molecules α, β, γ and bending angle θαβγ

実際の計算では，ある原子 α からカットオフ距離 (R + D) 内に存在する原子をリストアップし，カットオフ距離内の原子 β , γ , δ , … に対し，結合 α − β , β − γ , γ − δ ,… の二体間に働く力，及びそれぞれの結合から二つ，例えば結合 α − β と結合 α − γを選びその三体間によって発生する力をそれぞれ求め，これらの力のベクトルを重ね合わせることにより原子 α に働く力を求める．

Tersoﬀ ポテンシャルでは，T2(21) および T3(22) と呼ばれる 2 種類のパラメータセットが提案されている．表 2.1 にそれらの値を示す．T2 は，表面物性を優先させてフィッ
ティングしており，クラスターの形状を正確に表現することができるが (29) ，エネル
ギーが ab-initio の結果と異なることが報告されている (30) ．一方，T3 は，T2 では実験と一致しない弾性定数 C44 を修正し，結合角が 109.47◦ になるようにフィッティングしているため，バルクの記述には適しているが，結合角に依存した原子挙動を示すことが報告されている (29) ．本研究では T3 を用いた．

2.3 Tersoﬀ 型ポテンシャル

8

Table 2.1 Parameters for silicon.
Si − Si(T2)Si − Si(T3)3A [eV] 3.2647×10
1.8308×103
B [eV] 9.5373×101
4.7118×102−1
λ1 [Å ]3.23942.4799−1
λ2 [Å ]1.32581.7322−1
λ3 [Å ]1.32581.7322−1
β3.3675×10
1.1000×10−6n2.2956×101
7.8734×10−1c4.83811.0039×105d2.04171.6217×101
h0.0−5.9825×10−1
R [Å ]
3.0
2.85
D [Å ]0.2
0.15

原子に働く力は式 (2.2) で示すようにエネルギーの空間勾配によって表わされる．
[

]

∂E αβ
∂
r βα
αβ
αβ
αβ ∗ ∂
αβ
αβ
F =−
=
−
(f
(r
)f
(r
))
−
(b
)
(f
(r
)f
(r
))
c
R
c
A
∂r α
∂rαβ
∂rαβ
rαβ
∂
− fc (rαβ )fA (rαβ ) α (bαβ )∗
∂r
[
]
αβ
∂
r αβ
∂E
β
αβ
αβ
αβ ∗ ∂
αβ
αβ
F =−
=
−
(f
(r
)f
(r
))
−
(b
)
(f
(r
)f
(r
))
c
R
c
A
∂r β
∂rαβ
∂rαβ
rαβ
∂
− fc (rαβ )fA (rαβ ) β (bαβ )∗
∂r
αβ
∂E
∂
= − fc (rαβ )fA (rαβ ) γ (bαβ )∗(2.14)
Fγ = −γ∂r∂rα

ここで次の関係を用いた．

∂rαβ
∂rαβ
r βα ∂rαγ
∂rαγ
r γα
=
−
=
,
=
−
=
∂r α
∂r β
rαβ
∂r α
∂r γ
rαγ

(2.15)r αβ = r β − r α(2.16)F α ，F β の第一項は結合 α − β 間に働くの力である．これらの力は結合 α − β に関して独立に計算しておく．残りの項については結合 α − β と他の結合 (α − γ , α − δ ,…) との三体間を考慮して， (bαβ )∗ 及びその勾配を求めた上で計算する． (bαβ )∗ の勾2.3 Tersoﬀ 型ポテンシャル
配について以下に示す．

∂(bαβ )∗
1
=
∂r m
2

(

∂bαβ ∂bβα
+
∂r m
∂r m

9

)

(2.17)

1
∂bαβ
∂
1
=
{1 + β n (ζ αβ )n }− 2n −1{1 + β n (ζ αβ )n }
∂r m
2n
∂r m

(2.18)∂∂ζ αβ{1 + β n (ζ αβ )n } = nβ n (ζ αβ )n−1
∂r m
∂r m

(2.19)





∂  ∑
∂ζ αβ
=
fc (rαγ )g(θαβγ ) exp(λ33 (rαβ − rαγ )3 )
∂r m∂r m γ(̸=α , β)=∑∂fc (rαγ )g(θαβγ ) exp(λ33 (rαβ − rαγ )3 )∂r mγ(̸=α , β)(2.20)(m = α，β，γ)
ここで ζ αβγ = fc (rαγ )g(θαβγ ) exp(λ33 (rαβ − rαγ )3 ) とおくと，∂ζ αβγ ∂rαγ
∂ζ αβγ ∂(cos θ) ∂ζ αβγ ∂rαβ
∂ζ αβγ=++∂r α∂rαγ ∂r α∂(cos θ) ∂r α
∂rαβ ∂r α∂ζ αβγ r γα
∂ζ αβγ ∂(cos θ) ∂ζ αβγ r βα=++
∂rαγ rαγ∂(cos θ) ∂r α
∂rαβ rαβ()γα∂fc (rαγ )
3 αβ
αγ 2
αγ
3 αβ
αγ 3 r
−
3λ
(r
−
r
)
f
(r
)
g(θ)
exp(λ
(r
−
r
)
)
=
c
3
3
∂rαγ
rαγ
∂g(θ) ∂(cos θ)+fc (rαγ ) exp(λ33 (rαβ − rαγ )3 )∂(cos θ) ∂r α
r βα
+3λ33 (rαβ − rαγ )2 fc (rαγ )g(θ) exp(λ33 (rαβ − rαγ )3 ) αβr∂ζ αβγ
∂g(θ)
∂(cos
θ)
= fc (rαγ ) exp(λ33 (rαβ − rαγ )3 )β∂r∂(cos θ) ∂r β
r αβ
+3λ33 (rαβ − rαγ )2 fc (rαγ )g(θ) exp(λ33 (rαβ − rαγ )3 ) αβr)(αγ
∂ζ αβγ∂fc (rαγ )αγ 3 r
αγ 2
αγ
3 αβ
3 αβ
(r
−
r
)
)
(r
−
r
)
f
(r
)
g(θ)
exp(λ
=
−
3λ
c
3
3
∂r γ
∂rαγ
rαγ
∂g(θ) ∂(cos θ)+fc (rαγ ) exp(λ33 (rαβ − rαγ )3 )(2.21)∂(cos θ) ∂r γ(∂(cos θ)
cos θ
1
=
−
+
∂r α
rαβ
rαγ

)

(

cos θ
1
r βα
+
−
+
rαβ
rαγ
rαβ

)

r γα
rαγ

2.4 速度スケーリング法
(

10

)

cos θ r αβ
1 r αγ∂(cos θ)=−+∂r β
rαβ rαβrαβ rαγ()∂(cos θ)cos θ r αγ
1 r αβ= − αγ αγ + αγ αβ
∂r γ
r r
r r

(2.22)

また，∂g(θ)
2c2 (cos θ − h)= 2∂(cos θ)(d + (h − cos θ)2 )2(2.23)である．

2.4

速度スケーリング法

分子動力学法で温度制御する場合，もっとも簡単で直接的な方法として速度スケー
リング法がよく用いられる．熱統計力学より系の運動エネルギー K は次のように表さ
れる．

K=

N
3
1∑
mα (v α ·v α ) = N kB T
2 i=12

(2.24)

ここで，mα は原子 α の質量，v α は原子 α の速度，N は系の全原子数，kB はボルツマ
ン定数，T は系の温度である．式 (2.24) より，系の温度 T は原子速度を用いて，次のように求められる．
∑T =mα (v α ·v α )3N kB(2.25)

設定温度が TC ，式 (4.8) より求めたある時刻の温度が T のとき，速度スケーリング
法では，各原子の速度 v α を√

TC /T 倍し設定温度 TC に近づける．ベルレ法では，
2F∆r (t + ∆t) = r (t + ∆t) − r (t) = r (t) − r (t − ∆t) + (∆t)
α

√

を

α

α

α

α

α

(t)
(2.26)
m

TC /T ∆r α (t + ∆t) で置き換えることに相当する．平衡状態では，能勢の方法 (33)など外部との熱のやりとりをする変数を考慮した拡張系の分子動力学法によって得られるカノニカルアンサンブルに一致することが示されている．

2.5 高速化手法

2.5

11

高速化手法

領域分割による高速化
式 (2.5) からわかるように，N 個の原子からなる系では，Etot の評価に N ×(N −1) 回の原子対の計算が必要となる．一方，実際の結晶中では近接原子による遮蔽 (screening)
効果により第二近接距離程度より離れた原子はほとんど作用を及ぼさないことが知られている．このため，分子動力学計算では相互作用打ち切り (カットオフ) 半径 rc を導入し (図 2.2)，その半径内の原子からの寄与のみを考慮する．
しかしながら，相互作用する原子対の検索に N × (N − 1) 回の試行を要するため，系が大きくなるにつれ計算負荷が飛躍的に増加する．これを避けるために rc よりひとまわり大きい半径 rfc (図 2.2) 内の原子をメモリーに記憶し，rfc 内での原子対の探索とすることによりオーダー N の計算に近づける方法 (粒子登録法 (33) ) がこれまでよく用いられてきた．しかしながら，粒子登録法では rfc 半径より外の原子が rc 内に達すると力の評価が適切でなくなるので，一定のステップ毎に登録粒子の更新 (N × (N − 1)回の探査) を行わなければならない．このため，系がある程度の規模以上に大きくなると，粒子登録による高速化は登録更新の負荷により打ち消される．

rfc

rc

Fig.2.2 Schematic of bookkeeping method

領域分割法では，まず図 2.3 に模式的に示すようにシミュレートする系をカットオ
フ距離程度の格子状に分割する．ある原子に作用する力を評価する際には，その原子2.6 格子不安定性解析

12が属する領域（図 2.3 の着色部）と隣接領域内 (図 2.3 の斜線部) の原子からカットオフ
距離内の原子を探索する．原子が属する領域は，位置座標を領域ブロックの辺長 bx , by
で除した際の整数により判断できるので，領域分割そのものの計算負荷は小さい．領域分割法は，粒子登録法において登録更新の負荷が大きくなるような大規模な系の高速化に適している．

y
by
0

x
bx

Fig.2.3 Schematic of domain decomposition

2.6

格子不安定性解析

2.6.1

Born の安定基準

現在にいたるまで結晶の力学的安定性に関して様々な研究がなされているが，その基礎は格子力学にある (34) ．理解を容易にするため原著と表記を変えて説明するが，ひずみ εij の均一変形を受ける結晶の内部エネルギー U (εij ) は Taylor 展開により

1
U (εij )= U (0) + U ′ (0)εij + U (0)′′ εij εkl + · · ·2(2.27)と表すことができる．ここで，U ′ (0), U ′′ (0) はそれぞれひずみ εij による 1 次，2 次導関数を表す．結晶は εij = 0 において釣り合い状態にあるから，εij を変数とする 6 次元位相空間において，内部エネルギー曲面の勾配が 0，すなわち U ′ (0) = 0 でなければならない．ここで，ある状態 x における（無負荷とは限らない）応力 σij (x) と弾性係

2.6 格子不安定性解析

13

数 Cijkl (x) はそれぞれ

1 ∂U (x)V (x) ∂εij
1 ∂ 2 U (x)Cijkl (x) =V (x) ∂εij ∂εkl
σij (x) =(2.28)(2.29)と定義されることから，U ′ (0) = 0 の条件は x = 0 で応力が 0 であることに対応する．

V (x) は x における結晶の体積である．したがって，式 (2.27) において 3 次以降の高次項を無視して変形すると

1
U (εij ) − U (0) = V (0)Cijkl (0)εij εkl
2(2.30)

が導かれる．すなわち，弾性係数行列 Cijkl (0) が正定でなければ，∆U = U (εij ) − U (0)が負になり，x = 0 よりも低いエネルギー状態が存在することになる．

2.6.2

ブラベー格子を用いた安定性解析

上述の議論は理解しやすいが，式 (2.27) は無負荷状態における展開であることに注意が必要である．すなわち，外力下で釣り合いにある結晶の安定性も，単純に Cijkl (x)の正値性で評価できることにはならない．そこで，1970 年代に行われた格子不安定解析 (36)–(38) では，理想均一結晶の変形を，図 2.4 に示すように単位格子を記述する 6 つのパラメータ a1 ∼ a6 で代表させることで，外力下の結晶の安定性を議論している．すな
わち，内部エネルギー U を変数 ai について展開すると式 (2.27) と全く同じ形式になる
が，ひずみによる展開では 1 次微分の項が (2.28) 式の応力となり，外力下では 0 にならないのに対して，格子パラメータによる微分 U ′ (x) = ∂U (x)/∂ai は釣り合い状態にある任意の x で 0 となる．したがって，外力下の結晶の安定性は U ′′ (x) = ∂ 2 U (x)/∂ai ∂ajの正値性で議論される．この基準を用いて初めて，「応力−ひずみ経路の不安定分岐点」の存在が特定され，計算機能力が十分でない当時でも，[001] 単軸引張り下で現れる

α-Fe の bcc→bct 相変化などが解明されている (39) ．

2.6.3

弾性剛性係数による安定性解析

格子パラメータによる安定性解析では，原子の熱揺動を考慮することができず，有限温度下の結晶の安定性を議論することができない．そこで Wang, Yip, Li ら (40), (41)

2.6 格子不安定性解析

14Fig.2.4 Unit cell of fcc lattice.

は，応力と弾性係数によって表される弾性剛性係数 (42) の正値性によって結晶の安定性を評価する手法を提案した．現論文 (41) では有限ひずみ，温度下の結晶の自由エネ
ルギーを展開して導出しているが，単純に弾性剛性係数が，非線形弾性体の応力−ひずみ応答を表す物理量であることを考えれば，その正値性で結晶の安定性を評価しようと考えるのは自然な成りゆきといえる．具体的には，Voigt 対称性 (42) を持たせた次の弾性剛性係数 Bijkl の行列式の正負で安定性を評価する．

Bijkl = Cijkl + (σil δjk + σjl δik + σik δjl+ σjk δil − σij δkl − σkl δij )/2(2.31)

ここで δ はクロネッカーのデルタである．式 (2.29) の弾性係数は状態 x における内部エ
ネルギーにのみ依存するのに対し，弾性剛性係数は状態 x の外力（応力）の寄与が考慮される．応力，弾性係数は内部エネルギーの 1 次，2 次導関数で表されるが，Tersoﬀ型ポテンシャルは 3 体項があるため Cijkl は非常に複雑な形式となる (43) ．そこで本論
αを切断法によって求め，系の応力は
文では原子応力 σij

σij =

N
1 ∑
σα
N α ij(2.32)と求める．また前述のように，応力が 0 でない状態においてはひずみに対する応力変化が弾性剛性係数であるので，各原子を中心に周囲の原子配置を垂直方向に ±1 ％，せん断方向に ±0.5 ％の微小ひずみ摂動 (模式図 2.5) を与え，その時の応力勾配 ∆σij から数値的に評価した．

2.6 格子不安定性解析irc

Fig.2.5 Schematic view of local strain perturbation for evaluation of local
elastic stiﬀness.

15

第3章

[001] 引張シミュレーションならびに格子不安定解析

第一原理計算による「静力学的な」引張シミュレーションによって，単結晶シリコンの引張強度は εmax = 0.25, σmax = 18.73[GPa] と求められている (14) ．本章では，Tersoﬀ型ポテンシャル (21), (22) を用いたときの単結晶シリコンの [001] 方向引張における格子不安定条件を静力学解析で求めるとともに，分子動力学法で [001] 方向に単軸引張りシ
ミュレーションを行い，非弾性変形 (多数の原子の運動により可能な局所構造緩和) 発
生時の局所格子不安定について検討する．

3.1

静力学解析

3.1.1

解析条件

図 3.1 に示す Si ダイヤモンド構造の単位格子 (原子数 8) を用いて，その格子が無限に並んだ理想結晶についての解析を行う．格子定数は 0.543[nm] とした．[001] 方向にひずみ増分 ∆εzz =0.001(または，0.0001) を与えてアフィン変形 (z 方向の座標をスケー
リング) させ，そのときの応力を評価して応力−ひずみ曲線を求める．また，各ひずみ下で独立な 6 つのひずみ成分 (ε11 , ε22 , ε33 , ε12 , ε23 , ε31 ) 方向にひずみ摂動を与え，弾性剛性係数 Bij を求めた．なお，引張軸に垂直な方向については(i) x , y 方向にはスケーリングしない (ε´=0)16

3.1 静力学解析

17z [001]y [010]x [100]

Fig.3.1 Unit lattice of Si diamond structure.

(ii) 垂直応力が 0 となるように εxx = εyy の等方収縮条件で横方向にも座標スケーリン
グを行う (σ´=0)の 2 通りを考慮した．

3.1.2

応力−ひずみ関係

応力−ひずみ曲線を図 3.2 に示す．(i) の横ひずみ 0 での結果を黒色で，(ii) の等方収縮でのそれを赤色で表す．3 軸引張りに近い (i) では，応力−ひずみ曲線の勾配が σ´=0のそれより大きくなるが，fcc 金属の場合 (44) に比べて両者の差は小さい．(i) では，ひずみ εzz =0.378 でピーク応力 σzz =15.31[GPa] を示した後，ひずみ 0.397 からパルス状に応力がはねあがる不自然な曲線を示しているが，これは後述するように原子を固定した静力学で原子間隔が一方向に広げられた結果，別の結晶構造にスイッチしたためである．逆に σ´=0 の条件では εzz =0.595 において応力が不連続に急減している．この点では垂直応力が σ´=0 となる ε´が存在しなかったため，エネルギーが減少するように

ε´を変化させた結果，このように不連続な応力低下を生じた．この不連続を生じる直前の応力−ひずみ勾配は，正の値を有しており 0 ではない．εzz =0.595 の応力低下直前の最大応力は σzz =17.95[GPa] である．第一原理計算で求められている理想強度に比べ，引張ひずみは著しく大きくなったが最大引張応力は同程度となった．

3.1 静力学解析

Tensile stress σ zz [GPa]

20

15

10

5

ε' = 0
σ'=0
0

0.2
0.4
Tensile strain ε zz

0.6

Fig.3.2 Stress-strain curves under [001]tension (static analysis).

18

3.1 静力学解析

3.1.3

19

格子安定性

ε´=0 における detBij の変化を図 3.3 に示す．detBij はひずみの増加とともに単調に減少している．εzz ＞ 0.359 で detBij が負となっており，格子不安定点は εzz =0.36 と求められる．これは図 3.2 の応力ピークのひずみ 0.378 より小さいが，応力で見ると15.30[GPa] と両者に差はない．また，ひずみ 0.387 で再び detBij ＞ 0 となっていることから，別の安定な結晶構造にスイッチしたことがわかる．しかし，その後すぐに detBij
は急減し detBij ＜ 0 となっている．

σ´=0 での引張りにおける detBij の変化を図 3.4 に示す．0.51 までは detBij は滑らかに減少し，0.51 で初めて detBij ＜ 0 となった．ε´=0 のときに比べ，格子不安定となる
ひずみは最大応力を生じるそれ (0.595) より約 0.1 ほど小さくなったが，やはり第一原
理計算の結果と比べると大きな差がある．

x10 12Magnitude of 6x6 determinant, detBij , GPa 6

x10

12

0.02
4

2

0.01

0
0.32

0.34

0.36

0.38

0

0

0.1

0.2
0.3
Tensile strain ε zz

0.4

Fig.3.3 Relationship between tensile strain and crystal stability under [001]tension (ε´=0).

3.1 静力学解析

Magnitude of 6x6 determinant, detBij , GPa 6

x10
0.02

12

x10
1.5

0.015

7

1

0.5

0.01

0

0.005
0.45

0.48

0.51

0.54

0
0

0.2
0.4
Tensile strain ε zz

0.6

Fig.3.4 Relationship between tensile strain and crystal stability under [001]tension (σ´=0).

20

3.2 分子動力学シミュレーション

3.2

分子動力学シミュレーション

3.2.1

シミュレーション条件

21

図 3.1 の単位格子を，x, y, z 方向にそれぞれ 20 個ずつ並べた大きさ 10.86[nm]×10.86[nm]×10.86[nm] のシミュレーションセル (図 3.5) を解析対象とする．総原子数は 64000 個である．境界条件は (i) バルクシリコンを想定し，全方向周期境界条件としたもの，(ii)
自由表面からの欠陥生成を可能とするために，x, y 方向を自由境界条件としたもの，の

2 通りを考慮した．まず，それぞれの境界条件で 10000[fs] の初期緩和シミュレーショ
ンを行った．その後，[001] 方向に毎ステップひずみ増分 ∆εzz =5.0×10−6 を与えて引張りシミュレーションを行った．なお，(i) のバルクシリコンでは横方向の応力が 0 となるようにセル寸法をスケーリングしている．温度はいずれのシミュレーションにおy [010] 10.86 nm

10

z [001].86

nm10.86 nmいても 10[K] とし，速度スケーリングにより制御した．

x [100]
Fig.3.5 Simulation cell for tension.

3.2 分子動力学シミュレーション

3.2.2

22

シミュレーション結果と考察

全方向周期境界条件としたときの応力−ひずみ関係を図 3.6 に実線で示す．図には3.1.2 節の σ´=0 での静力学解析結果を破線で示している．欠陥がなく，また熱揺動の影響も小さいためひずみ εzz =0.411 まで両者は一致している．一方，先の格子不安定
ひずみ (εzz =0.51) より小さなひずみで非弾性変形を生じて応力が急減した．応力急減前後における原子配置を上から見たものを図 3.7 に，横から見たものを図 3.8 に示す．図中の原子はポテンシャルエネルギーの値に応じて着色している．応力急減前のひずみ εzz =0.411 では，原子のエネルギーはほぼ均一であるのに対して，応力急減後の εzz =0.412 では結晶が乱れた部分にエネルギーが高い，もしくは逆に低い原子が認められる．エネルギーが低い原子はエネルギーが高い原子に囲まれており，ループ状の欠陥生成 (転位に相当) とそれによるエネルギー緩和がなされたことを示唆している．

x, y 方向自由境界条件の場合の応力−ひずみ関係を図 3.9 に示す．先ほど同様に図3.1.2 節の σ´=0 での静力学解析結果を破線で示している．自由表面の存在により構造緩和の自由度が大きくなるため，静力学解析の応力−ひずみ曲線より低応力を示し，特に変形後期の高ひずみ域でそのずれが著しくなる．応力が急減したのはひずみ 0.388
で，(i) のバルクよりもさらに低ひずみ側で非弾性変形を生じた．応力急減前後における原子配置を図 3.10 ならびに図 3.11 に示す．先ほどと同様に，原子のエネルギーの値によって原子を着色し，上と横から見たものを表示している．自由表面の原子はダン
グリングボンドを有するため，結晶内部の原子に比べエネルギーが高く赤く着色されている．応力が低下するひずみ εzz =0.389 の図において，結晶表面の角部から欠陥が発生し，結晶内に転位が進行しはじめている様子が確認できる．

Tesnile stress σ zz

[GPa]

3.2 分子動力学シミュレーション

15

10

5

0

dynamics
statics

0.1

0.2
0.3
Tensile strain ε zz

0.4

Fig.3.6 Stress-strain curve under tension (bulk model, (i)).

0.5

23

3.2 分子動力学シミュレーション

[010]

[100]

(b) ε zz = 0.412(a) ε zz = 0.411

-4.1

Ε

-3.8[eV]Fig.3.7 Snapshots of atoms and energy distribution before and after the stressdrop (top view, bulk model).

[001]

[100]

(a) ε zz = 0.411-4.1(b) ε zz = 0.412

Ε-3.8[eV]Fig.3.8 Snapshots of atoms and energy distribution before and after the stressdrop (side view, bulk model).

24

3.2 分子動力学シミュレーション

Tensile stress σ zz

[GPa]

15

10

5

0

dynamics (wire)statics (bulk)

0.1

0.2
0.3
Tensile strain ε zz

0.4

Fig.3.9 Stress-strain curve under tension (ii).

0.5

25

3.2 分子動力学シミュレーション

[010]

[100]

(a) ε zz = 0.388-4.1(b) ε zz = 0.389

Ε-3.8[eV]Fig.3.10 Snapshots of atoms and energy distribution before and after the stressdrop (top view, wire model).

[001]

[100]

(a) ε zz = 0.388-4.1(b) ε zz = 0.389

Ε-3.8[eV]Fig.3.11 Snapshots of atoms and energy distribution before and after the stressdrop (side view, wire model).

26

3.2 分子動力学シミュレーション

3.2.3

27

局所格子不安定性

引張ひずみ下の「瞬間的な」原子配置データから，各原子を中心として周辺原子を独立な 6 つのひずみ成分 (ε11 , ε22 , ε33 , ε12 , ε23 , ε31 ) 方向にわずかに座標スケーリングを行い，そのときの原子応力変化から原子弾性剛性係数 Bijα (i, j = 1 ∼ 6) を数値的に求めた．その後，Bijα の 6×6 行列式の値 detBijα を計算し安定性判別を行った．先の図 3.7で示した，上から見たバルクの Si の原子配置を，detBijα ＜ 0 となった原子を赤色で着色して図 3.12 に示す．ひずみ 0.409 では不安定と判別された原子は存在しないが，ひずみ εzz =0.410 で初めて結晶内部に不安定原子が現れている．ひずみ 0.411 では，先の図 3.7 では結晶の乱れがなくポテンシャルエネルギーは一様であったが，局所格子不安定の基準でみると不安定原子が集合した領域が現れている (図 3.12(c))．εzz =0.410 の
図 (b) で最初に不安定原子が現れた領域が最も大きく成長しており，図 3.7 では実際にそこに欠陥が発生している．

x, y 方向を自由表面としたシリコンの原子配置を，同様に detBijα ＜ 0 の原子を赤色で着色して図 3.13 に示す．自由表面を持つ場合，引張前の平衡状態においても結晶表面に不安定原子が存在する (図 (a))．一方，引張ひずみを与えると表面の不安定原子が安定化し，結晶の角部 (エッジ) のみ不安定原子が存在する (図 (b))．εzz =0.387 において，エッジ部から結晶内部方向に不安定原子が増加している．ここで，図中で左上と1 ⃝1´) に多くの不安定原子が見られているのに対し，後で実際に欠右下のエッジ部 (⃝,2 ⃝2´) であり欠陥生成とともに不安定領域
陥が生成するのは右上−左下のエッジ部 (⃝,1 ならびに左下エッジ部 (⃝)
2 の原子構造
が結晶内部に進展している．左上エッジ部 (⃝)の詳細を，引張前 (εzz =0.0) と欠陥生成前の εzz =0.388 について図 3.14 および図 3.15にそれぞれ示す．図 3.1 の単位セルにおいて，着色した原子が実際に並べられる原子

1 のエッジは頂点に原子が存在するのに対し，⃝
2 のエッジは頂点には原子がなので，⃝1 の原子結合は引張による変化はなく「丸みをおびた」形状になっている．エッジ部 ⃝2 ではひずみの増加見られず，ダングリングボンドは変わらない．しかし，エッジ部 ⃝2 ⃝2´のエッジ部
によってエッジ部の原子がより結晶内部に移動している．そのため⃝,の方が不安定原子が少なかったものと考えられるが，こちらで欠陥を生じたのはエッジ頂点を「引き止める」原子が存在しなかったためと考えられる．

3.2 分子動力学シミュレーション(a) ε zz = 0.409(b) ε zz = 0.410(c) ε zz = 0.411(d) ε zz = 0.412

28

[010]

[100]

Fig.3.12 The distribution of the instability atoms under tension (bulk model, (i)).

3.2 分子動力学シミュレーション(a) ε zz = 0.00①

(b) ε zz = 0.350②'②(c) ε zz = 0.387

29①'(d) ε zz = 0.388

y [010]x [100](e) ε zz = 0.389
Fig.3.13 The distribution of the instability atoms under tension (wire model, (ii)).

3.2 分子動力学シミュレーション

y [010]x [100](a) ε zz = 0.0(b) ε zz = 0.388

1 in Fig.3.13 (d).
Fig.3.14 Detail of atomic configuration at edge ⃝y [010]x [100](a) ε zz = 0.0(b) ε zz = 0.388

2 in Fig.3.13 (d).
Fig.3.15 Detail of atomic configuration at edge ⃝

30

第4章
曲げシミュレーションならびに格子不安定解析

ナノスケールの単結晶シリコンビームの曲げ試験において，曲げ強度は温度に依存すること，高温で塑性変形が起こることなどが報告されている (1) ．本章では単結晶シ
リコンの無限平板およびナノワイヤの曲げシミュレーションを行い，寸法または温度による変形挙動の違いと，非弾性変形発生時の局所格子不安定性について検討する．

4.1

シミュレーション条件

4.1.1

無限平板モデル模式図 4.1 に示すように，スラブ状のシミュレーションセルに z 軸方向にのみ周期境界を適用して単結晶シリコンの無限平板をモデル化する．結晶方位は図に示したよ
うに x, y, z 軸をそれぞれ [100]，[010]，[001] としている．シミュレーションセルの寸法は表 4.1 に示すように Model I, II の 2 通りを考慮した．原子数はそれぞれ 115200 と14400 である．Model II は Model I の 1/2 のサイズとなっている．10000[fs] の初期緩和
1 ⃝
3 の表面原子に，変位 ∆d = 5.0 × 10−5 [nm]シミュレーションを行った後，図 4.1 の⃝∼を毎ステップ矢印方向へ増加することで三点曲げシミュレーションを行った．変位を与える原子はセルの z 方向幅によって変わり，Model I では 96，Model II では 24 である．温度は前章と同様 10[K] とし，速度スケーリングによって制御している．

31

4.1 シミュレーション条件

32

また，温度による影響を考慮するため，Model I のシミュレーションでは温度 300,

600[K] でも曲げシミュレーションを行った．

①

z [001]x [100]a1lhy [010]

a2

a3

w

②
③

Fig.4.1 Schematic model for bending simulation of infinite thin plate.

Table 4.1 Simulation model size.
Model I
Model IIl [nm]
65.16
32.58

h [nm] w [nm]
10.86
3.258
5.43
1.629

a1 [nm]
32.58
16.29

a2 [nm]5.43
2.715

a3 [nm]5.43
2.715

4.1 シミュレーション条件

4.1.2

33

ナノワイヤモデル

長さ 65.16[nm]，断面 10.86×10.86[nm2 ] の単結晶シリコンを解析対象とする (図 4.2)．
原子数は 384000 である．全方向自由境界の下，80000[fs] の初期緩和シミュレーション

1 ⃝
3 の表面原子 (それぞれ 400を温度 10[K] で行った．その後，先の無限平板と同様に⃝∼原子) に 5.0 × 10−5 [nm] の変位を矢印方向へ与え三点曲げシミュレーションを行った．
z [001]x [100]①
y [010]32.58[nm]65.16 [nm]10.86[nm]8.145[nm]②8.145[nm] 10.86[nm]③

Fig.4.2 Model for bending simulation of nano wire.

4.2 シミュレーション結果と考察

4.2

シミュレーション結果と考察

4.2.1

無限平板の変形挙動

34

1 の表面原子に作用する
Model I(T=10[K]) のシミュレーションにおいて，図 4.1 の⃝
2 ⃝
3 の表面原子位置を基準とした⃝
1 の相対移動距離から変位を算出
y 方向の力と，⃝,して力−変位曲線を求めた (図 4.3)．変位速度が大きいため，応力波によるのこ歯状の応答を示しているが，全体的には押し込み量の増加とともに押し込み力が増加している．押し込み変位 d=10[nm] 過ぎでは押し込み力の振動が減衰しているようにみえるが，d=23.8[nm] から再び振幅が大きくなり，最終的に d=39.0[nm] で平板がぜい性的に割れて押し込み力が急減した．またわずかではあるが，詳細に見ると変位 d=23.5[nm]の点 (図 4.3(d) 近傍) で反力が一時的に減少している．図 4.3 で (a)∼(h) で示した押し込

1 の点近傍を拡大して図 4.4 に示す．図ではポテンシャル
み変位における原子配置を，⃝
エネルギーの値に応じて原子を着色している．前章と同様，変位 d=0.0[nm](a) の図ではダングリングボンドを有する自由表面原子のみ赤く着色されている．曲げモーメン
トの増加によって，曲げ応力が大きくなる部分にエネルギーの高い原子が増加している(図 4.4(b))．変位 d=20.0[nm] の図 4.4(c) では，変位を与えた表面原子が結晶内に押し込まれ，その下に高いエネルギーを持つ原子 (黄色) が認められる (図 4.4(c), (d))．(c)から (d) の一時的な力の減少は表面原子の押し込みによるステップ形成のためと考えられる．その後，変位を与えた原子の両端から結晶内部へ転位が発生している (図 4.4(e))．
このとき，引張側では欠陥はまだ見られない．力が急減しはじめる変位 d=37.2[nm] から，それまで大きな変化がなかった引張応力側の表面からぜい性的に破断している (図4.4(h))．

4.2 シミュレーション結果と考察

0.2(f), (g)(c)

0.15Force fy [ µ N]

35

0.1

(d)

0.05

(e)

(b)

(a)
(h)

0
0

10

20

30

Displacement d [nm]Fig.4.3 Force-displacement curve (Model I, 10[K]).

40

4.2 シミュレーション結果と考察(a) d = 0.0 [nm](e) d = 30.0[nm](b) d = 10.0 [nm](f) d = 37.1[nm](c) d = 20.0 [nm](g) d = 37.2 [nm](d) d = 25.0 [nm](h) d = 39.0 [nm]

-4.6

Ε

-3.8[eV]Fig.4.4 Snapshots of atoms and energy distribution (Model I, 10[K]).

36

4.2 シミュレーション結果と考察

37

寸法効果
寸法の小さい Model II の力−変位曲線を図 4.5 に示す．図中には Model I の結果も細線であわせて示している．原子数が少ないので応力波振幅 (ノイズ) は小さい．また変位制御で押し込みを行っているので，反力の総和は Model I に比べ小さくなっている．Model I と同様，力の振幅が小さくなった後，記号 (d) で示したように一時的な力の減少とともに再び振動している．(f), (g) で示した変位 d=20.1[nm] で力が急減する
が，Model I と異なり力は 0 にならない．
図 4.5(a)∼(h) で示した点における原子配置を図 4.6 に示す．力が一時的に減少した
点 (d) では，やはり変位を与えた表面原子が結晶内に押し込まれステップを形成している．全体的な力の低下 (d=20.1[nm] 以降) はやはり引張側で非弾性変形を生じたときであるが，寸法の小さい Model II では Model I のときのようにぜい性的に破断することなく表面側に乱れが導入されるだけで変形は停止した．これは系が小さいために，
たくわえられた弾性ひずみエネルギーが小さいためと考えられる．

0.2

Force fy [ µ N]

0.15

0.1

(c)

(e)

0.05

(a)

(b)

0
0

10

(d)
20

Model I
Model II(f),(g)(h)30

40Dispalacement d [nm]Fig.4.5 Force-displacement curves under bending(Model II).

4.2 シミュレーション結果と考察(a) d = 0.0 [nm](e) d = 16.0 [nm](b) d = 5.0 [nm](f) d = 20.1 [nm](c) d = 10.0 [nm](g) d = 20.2 [nm](d) d = 13.7 [nm](h) d = 20.6 [nm]

-4.6

Ε

-3.8[eV]Fig.4.6 Snapshots of atoms and energy distribution(Model II).

38

4.2 シミュレーション結果と考察

39

温度の影響

Model I の 300, 600[K] における力−変位関係をそれぞれ図 4.7，図 4.8 に示す．10[K]の場合と比べて揺らぎが大きいが，これは応力波だけでなく原子の熱揺動の寄与も含まれる．温度が高くなるにつれ，押し込み後期の反力が小さくなり，また力が急減し0 となる変位も 300[K] で 27.1[nm]，600[K] で 25.6[nm] と小さくなっている．
図 4.7，図 4.8 の (a)∼(h) における原子配置をそれぞれ図 4.9 ならびに図 4.10 に示した．いずれも最終的な応力急減は，引張側表面からのぜい性的な破断によりもたらさ
れているが，詳細にみると力が低下する直前の (e) 点で引張を受ける側の表面の○で囲んだ箇所で原子配列に乱れが生じはじめている．(f) の点でもまだ押し込み反力が急減していないが，上下表面に段差を生じており，わずかではあるが「塑性変形」している．実際，図 4.7，図 4.8 の力−変位曲線には，反力が増加せずに変位が増加しているプラトー領域が存在する．

0.2

0.15

Force fy [ µ N]

(f)
(b)
0.1

(e)

0.05

(c)

(d)
(g)

(a)

(h)

0
0

10

20

30

Displacement d [nm]Fig.4.7 Force-displacement curve under bending(300[K]).

40

4.2 シミュレーション結果と考察

0.2

Force fy [ µ N]

0.15

(c)
(f)

(b)
0.1

0.05

(d) (e)
(g)

(a)

(h)

0
0

10

20

30

Displacement d [nm]Fig.4.8 Force-displacement curve under bending(600[K]).

4040

4.2 シミュレーション結果と考察(a) d = 0.0 [nm](e) d = 25.2 [nm](b) d = 10.0 [nm](f) d = 27.8 [nm](c) d = 15.0 [nm](g) d = 28.0 [nm](d) d = 20.0 [nm](h) d = 28.5 [nm]

-4.6

Ε

-3.8[eV]Fig.4.9 Snapshots of atoms and energy distribution (300[K]).

41

4.2 シミュレーション結果と考察(a) d = 0.0 [nm](e) d = 23.8 [nm](b) d = 10.0 [nm](f) d = 24.5 [nm](c) d = 15.0 [nm](g) d = 25.0 [nm](d) d = 20.0 [nm](h) d = 25.8 [nm]

-4.6

Ε

-3.8[eV]Fig.4.10 Snapshots of atoms and energy distribution(600[K]).

42

4.2 シミュレーション結果と考察

4.2.2

43

ナノワイヤの変形挙動

ナノワイヤを押し込んだ時の力−変位関係を図 4.11 に示す．押し込み開始直後の振動が，先の平板の曲げに比べて大きくなっている．押し込み量 d=15.5[nm] で最大反力 fy =0.9[µN] を示してその後急減した．図 4.11 で示した (a)∼(g) 各点における原子配置を図 4.12 に示す．図 4.12 左上の模式図のように，最大曲げモーメントを生じる部分(赤色破線で囲まれた箇所) を斜め下方向から見ている．自由表面があっても，ワイヤ全体を横断するような転位の運動は認められず，平板と同様引張側表面からのぜい性的な破断を生じている．

1(d),(e)0.8

Force fy [ µ N]

(c)
0.6

(f)
0.4

(g)
0.2

(b)

(a)
0

10

20

Displacement d [nm]Fig.4.11 Force-displacement curves under bending(nano wire).

4.2 シミュレーション結果と考察(d) d = 15.5 [nm]Schematic view(a) d = 0.0 [nm](e) d = 15.6 [nm](b) d = 5.0 [nm](f) d = 16.0 [nm](c) d = 10.0 [nm](g) d = 18.7 [nm]

-4.6

Ε

-3.8[eV]Fig.4.12 Snapshots atoms and energy distribution of bending(nano wire).

44

4.2 シミュレーション結果と考察

4.2.3

45

局所格子不安定性

無限平板

Model I, 10[K] のシミュレーションにおける「瞬間的な」原子配置を用いて局所格
子不安定解析を行い，detBijα の正負に応じて着色したものを図 4.13 に示した．図では

2 ，⃝3 の制御部分でも結
detBijα ＜ 0 となった原子を赤色で着色している．3 点曲げの⃝1 ⃝
3 の押し込み部分近傍に不安定原子が確認される (図晶内部への押し込みを生じ，⃝∼
3 の支持部近傍に垂直方向に並んだ不安定原子が認められる．さ4.13(b))．図 (c) では⃝
2 ，⃝3 の支持部の周囲が全体的に detBijα ＜ 0 となっている．逆に曲げらに図 (d) では⃝
中心部には detBijα ＜ 0 になった原子は少なく，圧縮側の欠陥生成部にいくつか認められるだけである．図 (c) の直線状の不安定原子から推測するに，板のせん断および回転に対して不安定となったものと考える．板が大きく U 字状になった (f) 点では，これらの不安定原子の領域は変位を与えた点の中心に遷移している．割れを生じる引張側の表面には確かに不安定原子が現れているが (図 4.13(f))，全ての弾性剛性係数成分を考慮した detBijα ＜ 0 の条件では，先述の全体的な不安定原子と，最終的に破断を生じるそれを区別することはできない．

Model II の原子配置についても同様に，detBijα ＜ 0 となった原子を赤色で着色して
2 ，⃝
3 の支持
図 4.14 に示した．基本的には Model I と同じような変化を示しており，⃝部分に不安定原子が多数発生しそれが中央部分に移動している．Model I との違いは，

1 の引張表面側にも不安定原子が多く見られることであり，力急減時には実際にこの⃝部分に非弾性変形を生じている．

Model I で温度を上げた時の結果を図 4.15 ならびに図 4.16 に示す．
「瞬間的な」原子配置で安定・不安定を判別しているため，10[K] の場合に比べ，力学的条件だけでなく

2 ，⃝
3 近傍の不安定原子「確率的な」ゆらぎを含んだ判定となっている．先述の支持部⃝
群は，300[K] ではその傾向が見られるが，600[K] ではほとんど見られなくなっている(図 4.15(d)∼(h), 図 4.16(e)∼(h))．また，300[K] の場合も，これらの不安定原子は中心部に移動することはない．その理由としては，10[K] のときほど大きく U 字状に曲げられる前に破断したことがあげられる．また破断部には不安定原子が見られるが，温4.2 シミュレーション結果と考察

46

度が高くなるとその傾向ははっきりせず，瞬間的な detBijα の正負でその発生を予測するのは難しい．
ナノワイヤ
図 4.17 に detBijα ＜ 0 の不安定原子を赤色で着色した原子配置を示す．引張前の図(a) では表面原子が赤く着色されているが，押し込みを行うと引張を受ける表面の不安
1 ⃝
3 付近以外では安
定原子は安定化する (図 4.17(b))．さらに側面でも変位を与えた⃝∼
2 ，⃝3 近傍では内部の原子も不安定化している (図 4.17(c))．これまでと同様，支持部⃝
1 の押し込み部では側面に不安定原子が残留している．
定となる．無限平板と異なり，⃝a にはエッジ部にわずかに不安定原
欠陥生成直前の図 4.17(d) において，割れて生じる⃝子が存在するだけである．
以上のように，Bijα の 6×6 の全ての成分を考慮した条件では，最終的な不安定破壊の予測につながらない．そこで，弾性剛性係数成分の主応力部分 (6×6 行列の中の 3×3 主小行列，対角化したときの σ1 , σ2 , σ3 を含む成分) が負となった原子を無限平板 Model

I，T=10[K] の場合について調べ，赤色で着色して図 4.18 に示した．主応力成分の判定
2 ，⃝3 の支持部の不安定原子群は見られず，これらの不安定がせん断または回では，⃝転に対するものであることがわかる．表面原子は依然として不安定と判断された原子が存在するが，最終破断を生じる変位 d=37.1[nm] の図 (f) では集団的に不安定となった領域を生じている (図 4.18(e))．他のモデルでも同様な傾向が見られた．
理想結晶に対する静力学解析では，静水圧応力が結晶の安定限界を上げることが示されている．したがって，例えば図 4.13(f) の同じ黄色で着色した「安定な」原子にお2 ，⃝3 の近傍と⃝1 の下の部分ではその力学状態は大きく異なる．瞬間
いても，支持部⃝の安定・不安定だけでなく，安定 → 不安定となった時の「カタストロフ」さを考慮す
ればより正確な欠陥・破壊発生予測につながるものと思われるが，現時点ではその方法は不明である．

4.2 シミュレーション結果と考察①③②(a) d = 10.0 [nm](e) d = 33.0 [nm](b) d = 20.0 [nm](f) d = 37.1 [nm](c) d = 25.0 [nm](g) d = 37.2 [nm](d) d = 30.0 [nm]detBijα > 0(h) d = 39.0 [nm]αijdetB < 0

Fig.4.13 The distribution of the instability atoms in the infinite thin plate
under bending (Model I, T=10[K]).

47

4.2 シミュレーション結果と考察(a) d = 5.0 [nm](e) d = 20.1 [nm](b) d = 10.0 [nm](f) d = 20.2 [nm](c) d = 15.0 [nm](g) d = 20.6 [nm]detBijα > 0
αdetBij < 0(d) d = 16.4 [nm]Fig.4.14 The distribution of the instability atoms in the infinite thin plate
under bending (Model II, T=10[K]).

48

4.2 シミュレーション結果と考察(a) d = 10.0 [nm](e) d = 26.0 [nm](b) d = 15.0 [nm](f) d = 27.0 [nm](c) d = 20.0 [nm](g) d = 28.0 [nm]α(d) d = 25.0 [nm]detBij > 0(h) d = 28.5 [nm]αijdetB < 0

Fig.4.15 The distribution of the instability atoms in the infinite thin plate
under bending (Model I, T=300[K]).

49

4.2 シミュレーション結果と考察(a) d = 10.0 [nm](e) d = 22.0 [nm](b) d = 15.0 [nm](f) d = 23.0 [nm](c) d = 20.0 [nm](g) d = 24.0 [nm]α(d) d = 21.0 [nm]detBij > 0(h) d = 25.0 [nm]αijdetB < 0

Fig.4.16 The distribution of the instability atoms in the infinite thin plate
under bending (Model I, T=600[K]).

50

4.2 シミュレーション結果と考察

51(c) d = 10.0 [nm]Schematic viewa(d) d = 15.5 [nm](a) d = 0.0 [nm]bα(b) d = 5.0 [nm]detBij > 0α(e) d = 15.6 [nm]detBij < 0

Fig.4.17 The distribution of the instability atoms in the bending transformation (nano wire).

4.2 シミュレーション結果と考察(a) d = 10.0 [nm](e) d = 33.0 [nm](b) d = 20.0 [nm](f) d = 37.1 [nm](c) d = 25.0 [nm](g) d = 37.2 [nm]α(d) d = 30.0 [nm]detBij > 0(h) d = 39.0 [nm]αijdetB < 0

Fig.4.18 The distribution of the instability atoms in the bending transformation (Model I, T=10[K]).

52

第5章
ナノカンチレバーの曲げシミュレーションならびに格子不安定解析

シリコンカンチレバーの押し込み試験では，結晶方位によって変形挙動が異なることが報告されている (6) ．本章では，単結晶で異なる結晶方位としたカンチレバーなら
びに bamboo 状に粒界を有するカンチレバーの押し込みシミュレーションを行い，その変形挙動と局所格子不安定性について検討する．

5.1

シミュレーション条件

解析モデルを模式図 5.1 に示す．カンチレバー部は長さは 32.58[nm]，断面 5.43×5.43[nm2 ] であり，土台は 16.29[nm]×21.72[nm]×21.72[nm] である．カンチレバー及び土台
1 ダイヤモンド構造の Si 単結晶から切り出す (カンチレバー部の x, y, z の結晶方位は⃝2 カンチレを ([100], [010], [001])，([110], [1̄10], [001])，([111], [1̄10], [1̄1̄2]) の 3 通り)，⃝
バー部で軸方向の結晶方位を長さ 5.43[nm] ごとに先端から [100]，[110] と交互に変え

1 の [100] に対応) の 4 通りを考慮した．温度は 10[K] とたバンブー構造とする (土台は⃝して，速度スケーリングにより制御した．まず，全方向自由境界条件の下で 80000[fs]の初期緩和シミュレーションを行った．その後，カンチレバーを支える土台の表面原

a にある表面原子に矢印方子を拘束し，カンチレバー先端から 1.5[nm] の位置 (図 5.1⃝)
向 (y 軸方向) に変位 1.0 × 10−4 [nm] を毎ステップ与えて押し込みシミュレーションを行った．各モデルの原子数ならびに変位を与えた原子の数を表 5.1 に示す．

53

5.1 シミュレーション条件

16.29 [nm]32.58 [nm]a5.42[nm]
5.42[nm]
5.42[nm]

5.42[nm]
21

[n
.72

z

m]

x
y

Fig.5.1 Simulation model of nano cantilever.

Table 5.1 Simulation condition.
silicon
[100]
[110]
[111]
bambooNo. of atoms
No. ofunder
all atoms
displacement control
1343446012414425143968431238788054

5.2 シミュレーション結果と考察

5.2

シミュレーション結果と考察

5.2.1

押し込み荷重−変位関係

55

変位を制御した表面原子に作用する y 方向の力を合計して，カンチレバーにかかる力

1 の結晶方位 [100]，[110]，[111] の単結晶シリコン，⃝
2 のバンブー構造として評価した．⃝
へそれぞれ下向きに変位を与えて押し込みを行ったときの力−変位関係を図 5.2(a)∼(d)に示す．力−変位関係では単結晶とバンブー構造に大きな差は無い．押し込んだ表面原子に与えた変位速度がやや大きいために表面原子が結晶内に押し込まれ，その反力によって押し込み初期に大きな応答を示した．力の振幅が大きいため見にくいが，押し込み量増加とともに押し込み力は徐々に増加している．しかし，力の急減などの大きな変化は見られなかった．各変位のときのスナップショットをカンチレバー根元付
近を拡大して結晶方位 [100] の単結晶を図 5.3，[110] を図 5.4，[111] を図 5.5，そしてバ
ンブー構造を図 5.6 に示した．これまでと同様，ポテンシャルエネルギーの値に応じて原子を着色している．単結晶シリコンでは，結晶方位によって表面およびエッジ部のエネルギーに違いがあるものの変形挙動に大きな差は見られなかった．バンブー構造
では粒界部分にエネルギーが高い原子が存在し，そのためカンチレバーつけ根部分の変形が単結晶に比べてひずみ (ポテンシャル) エネルギーが大きいが，それ以外に大きな違いは見られない．押し込み量が増加するにつれ，どのカンチレバーでも引張または圧縮を受ける側の結晶内部でエネルギーが高くなっているが，単軸引張または曲げシミュレーションで欠陥生成したときに見られた高いエネルギー集中は見られず，結晶の原子配列に乱れは生じなかった．d=20.0[nm] のときのカンチレバー上側の最大引
張応力，下側の最大圧縮応力を表 5.2 に示す．単結晶では，最大応力はいずれもカンチ
レバーの根元の，圧縮を受ける下側のエッジ部で生じている．バンブー構造では，先述のようにカンチレバーの根元部分に粒界があり，この粒界原子の引張側部分の最大引張応力が，下側の最大圧縮応力より大きくなった．前章のナノワイヤのシミュレー
α=18.4[GPa]
ションにおいて，最終的なぜい性破壊を生じるときの最大引張応力は σxxであったのに対し，本章で観察した原子応力はいずれもこれより小さい．カンチレバー

5.2 シミュレーション結果と考察

56

a までの距離が前章図 4.2⃝∼
1 ⃝
3 の距離より大きく，またナの根元から，変位を与えた⃝
ノワイヤに比べ幅，厚さが 1/2 であったため曲げによって蓄えられる弾性ひずみエネ
ルギーが小さく，d=20.0[nm] の曲げでも欠陥生成に到らなかったものと考える．

Table 5.2 Number of atoms

[100]
[110]
[111]
bambooMaximum tensilestress [GPa]
11.42
13.81
13.71
17.02

Maxmum compression
stress [GPa]-13.60-16.29-17.18-15.37

5.2 シミュレーション結果と考察

0.04

[µ N]

[ µ N]

0.04

57

Force

fy
Force

0

-0.02

0.02

fy

0.02

0

10

-0.02

20

Displacement

0

0

[nm]

d

20

Displacement

[nm]

d

(b) Single [110](a) Single [100][µ N][ µ N]

0.04

0.04

0.02

Force

fy

fy

0.02

Force

10

0

-0.02

0

10

-0.02

20

Displacementd(c) Single [111]

0

[nm]

0

10

Displacement(d) Bamboo

Fig.5.2 Force-displacement curves.

20

d

[nm]

5.2 シミュレーション結果と考察(a) d = 0.0 [nm](d) d = 15.0 [nm](b) d = 5.0 [nm](e) d = 20.0 [nm]

-4.6

Ε

-3.8[eV](c) d = 10.0 [nm]1
Fig.5.3 Snapshots of atoms and energy distribution (⃝single[100]).

58

5.2 シミュレーション結果と考察(a) d = 0.0 [nm](d) d = 15.0 [nm](b) d = 5.0 [nm](e) d = 20.0 [nm]

-4.6

Ε

-3.8[eV](c) d = 10.0 [nm]1
Fig.5.4 Snapshots of atoms and energy distribution (⃝single[110]).

59

5.2 シミュレーション結果と考察(a) d = 0.0 [nm](d) d = 15.0 [nm](b) d = 5.0 [nm](e) d = 20.0 [nm]

-4.6

Ε

-3.8[eV](c) d = 10.0 [nm]1
Fig.5.5 Snapshots of atoms and energy distribution (⃝single[111]).

60

5.2 シミュレーション結果と考察(a) d = 0.0 [nm](d) d = 15.0 [nm](b) d = 5.0 [nm](e) d = 20.0 [nm]

-4.6

Ε

-3.8[eV](c) d = 10.0 [nm]2
Fig.5.6 Snapshots of atoms and energy distribution (⃝bamboo).

61

5.2 シミュレーション結果と考察

5.2.2

62

局所格子不安定性

変位 d=20.0[nm] のときのカンチレバーのスナップショットを，格子不安定性を考慮して図 5.7，図 5.8 にまとめて示す．図 5.7 は引張表面側を，図 5.8 は圧縮表面側を見たものである．Bijα の主応力成分のみ考慮し，3×3 行列式の値が負となった原子を赤色で着色している．土台部分を見ると，3×3 の主小行列式を用いても {001} 表面は不安定と判定されるのに対し，高指数面である {1̄1̄2} 表面の方が不安定原子が少ない．また，
図 5.7(b) の [110] モデルは図 5.7(a) と同じ {001} 表面であるにもかかわらず，行列式が負になった原子が少ない．6×6 行列式での判定では，[110] モデルの土台側面全体にも不安定原子が見られた．これは，z 方向の結晶方位が同じ [001] であっても x, y 方向の方位が異なるとき，ε14 , ε15 , ε16 などの成分が無視できないことを意味する．カンチ

1
レバー部に着目すると，図 5.7，図 5.8(a) の⃝[100]では，カンチレバーの上側または下側の表面原子は一度全体に安定化した後，引張側ではカンチレバー中心領域，圧縮
側ではカンチレバーつけ根近傍に再び不安定となる原子が現れる．しかしながら，前章の曲げシミュレーションの不安定破断前に見られたような，結晶内部での不安定領

1
域は見られなかった．⃝[110]，[111] では，カンチレバーの側面において不安定原子が現れたが，[100] で見られた上側または下側で不安定原子は現れなかった (図 5.7(b)(c),
図 5.8(b),(c))．バンブー構造では，カンチレバーの引張を受ける側で [100] で見られた領域で不安定原子が見られた (図 5.7(d))．また，結晶方位が異なる粒界では表面と粒界の「交線」で不安定原子の集団が見られるが，中心の粒界面では現れなかった．
カンチレバーの曲げにおいては，変位 d=20.0[nm] を与えても表面の原子で安定・不安定の変化は見られたが，結晶内では不安定原子は現れなかった．また，先述のように，3×3 行列の成分だけの判定では座標軸に依存するので，本来は 6×6 行列を対角化しその固有値の正値性を議論することが望ましい．これらが今後検討されることを期待する．

5.2 シミュレーション結果と考察(a) single [100](c) single [111](b) single [110]detBijα > 0(d) bamboo

αdetBij < 0
Fig.5.7 The distribution of the instability atoms (top view, d=20.0[nm]).

63

5.2 シミュレーション結果と考察(b) single [110](a) single [100](c) single [111]detBijα > 0(d) bamboo

αijdetB < 0
Fig.5.8 The distribution of the instability atoms (bottom view, d=20.0[nm]).

64

第6章
結論

本研究は，分子動力学法を用いて様々な形状のシリコンについて変形シミュレーショ
ンを行うとともに，非弾性変形開始時の力学状態を局所格子不安定性に基づいて検討した．以下に，得られた結果を総括する．
第 2 章では，解析手法の概要を述べた．はじめに，分子動力学法の基礎方程式を示
し，原子間相互作用の評価に用いた Tersoﬀ 型ポテンシャルの概要を説明した．次に，分子動力学計算における温度制御方法と大規模分子動力学計算に用いられる領域分割による高速化について述べた．さらに，各原子の安定性を議論する局所格子不安定性
解析について概説した．
第 3 章では，用いた Tersoﬀ 型ポテンシャル (21), (22) について，理想完全結晶状態で [001]方向に単軸引張を与えたときの弾性安定限界を静力学解析で求めるとともに (局所じゃな
い従来の格子不安定解析)，分子動力学法で [001] 方向に単軸引張りシミュレーションを行い，多数原子の存在によって可能な非弾性変形発生を観察するとともにその時の局所不安定性について検討した．まず静力学解析では，x, y 方向にはスケーリングしない場合(ε´=0) では応力−ひずみ曲線はひずみ εzz =0.378 でピーク応力 σzz =15.31[GPa] を示すが，格子不安定となるひずみはそれより小さく εzz =0.36 である．ただし，そのときの応力は15.30[GPa] であまり差はない．垂直応力 0 となるように横方向等方収縮条件で座標ス
ケーリングを行った場合 (σ´=0) では，ひずみ εzz =0.595 で他の構造にシフトするため応力−ひずみ曲線が不連続に急減するが，それより小さいひずみ 0.51 で detBij ＜ 0 となった．ただし，ε´=0，σ´=0 いずれも第一原理計算の結果 (εmax = 0.25, σmax = 18.73[GPa])65

66と比べると格子不安定となるひずみは Tersoﬀ ポテンシャルでは過大評価されることがわかった．ただし，応力で見ると大きな差はない．分子動力学法によるバルクシリ
コンの引張シミュレーション (σ´=0) ではひずみ εzz =0.411 で結晶内に乱れを生じ，応力が急減した．引張ひずみ下の「瞬間的な」原子配置データから detBijα ＜ 0 となった
「不安定」原子を調べると，応力急減直前のひずみ 0.410 で初めて結晶内に不安定原子が現れ，そこを中心として不安定原子領域が拡大してループ状の欠陥が生成していた．
また，x, y 方向を自由境界条件としたナノワイヤの引張シミュレーションでは，引張ひずみ 0.389 でエッジ部から転位が発生して応力が急減したが，やはり応力急減前にエッジ部に不安定原子が現れていた．ただし，実際に欠陥が発生したのは 4 つあるワイヤエッジ部で「後から」不安定原子が発生した 2 つのエッジで，先に不安定原子を生じた他の 2 つのエッジは多くの不安定原子を生じているのに欠陥生成に到らなかった．この違いは，エッジ部に Si 原子列が存在するか否かによることを明らかにした．
第 4 章では，無限平板およびナノワイヤの曲げシミュレーションを行い，寸法または温度による変形挙動の違いと，非弾性変形発生時の局所格子不安定性について検討した．無限平板の寸法が大きいモデルでは，温度によらず引張を受ける側の表面にぜい性的な破断を生じて力が急減した．小さいモデルでは蓄えられる弾性ひずみエネル
ギーが小さく，引張側表面に非弾性変形を生じて力が低下したが，力が 0 になる破断は生じなかった．ナノワイヤの曲げでも，10[K] の条件ではワイヤ全体を横断するよ
うな転位は見られず，平板同様にぜい性的な破断を生じた．これらのシミュレーショ
ン中の原子配置から，各原子の Bijα の正値性を評価した結果，全ての原子弾性剛性係
数成分を考慮した 6×6 行列の detBijα ＜ 0 の条件では，変位を与えた両端の支持部近傍で板の回転およびせん断による大域的な不安定原子群が現れ，最終的に破断を生じる部分と区別することができなかった．そこで 6×6 行列の主小行列である 3×3 行列(ε11 , ε22 , ε33 , ε12 , ε23 , ε31 に関する部分) の正値性を見たところ，上述のように変形を生じない部分での不安定原子がほとんど見られなくなり，最終的に破断を生じる部分が破断直前に detBijα (3×3) が負となることがわかった．
第 5 章では，軸方向の結晶方位が異なる単結晶のカンチレバーならびに bamboo 状に結晶を有するカンチレバーの押し込みシミュレーションを行い，その変形挙動と局67
所格子不安定性について検討した．長さ 48.87[nm] のカンチレバー先端を d=20.0[nm]まで押し込んだが，いずれもカンチレバー反力が急減することはなく，また結晶方位
または結晶構造による違いも見られなかった．原子配置を見ても初期構造からの大きな乱れはなく，非弾性変形を生じていなかった．このため前章同様，原子弾性剛性係数の 3×3 行列の正値性を見ても，結晶表面には不安定原子が見られたが，結晶内部に不安定原子が現れることはなかった．なお，3×3 行列による判定では応力評価時の座標
軸のとり方に依存するため，6×6 行列の固有値解析を行う方が望ましいことも示した．
局所格子不安定性解析では，瞬間的な原子配置において各原子を「安定」「不安定」を議論する．したがって，無負荷でのつりあいと，臨界外力下でのつりあいいずれにおいても原子は同じ「安定」と判断される．今後は，瞬間の安定・不安定だけでなく，安定 → 不安定となった時の「カタストロフ」さを考慮したより正確な欠陥・破壊発生予測手法の検討が望まれる．

参考文献(1) Takahiro Numazu, Yoshitada Isono, Member IEEE, and Takeshi Tanaka, Plastic
Deformation of Nanometric Single Crystal Silicon Wire in AFM Bending Test
at Intermediate Temperatures, Journal of Microelectromechanical Systems 11(2002) 125.
(2) Hsien-Kuang Liu, C. H. Pan, Pang-Ping Liu, Dimension eﬀect on mechanicalbehavior of silicon micro-cantilever beams, Measurement 41 (2008) 885.
(3) Xueping Li, Takashi Kasai, Shigeki Nakao, Taeko Ando, Mitsuhiro Shikida,
Kazuo Sato, Hiroshi Tanaka, Anisotropy in fracture of single crystal silicon film
characterized underuniaxial tensile condiction, Sensors and Actuators A 117(2005) 143.
(4) 土屋智由, 薄膜の引張強度評価, 豊田中央研究所 R&D レビュー 34 (1999).
(5) Yuhang Jing , Qingyuan Meng and Wei Zhao, Molecular dynamics simulationsof the tensile and melting behaviours of silicon nanowires, Physica E 41 (2009)685.
(6) Yuhang Jing , Qingyuan Meng , Yufei gao, Molecular dynamics simulations of
the buckling behavior of silicon nanowires under uniaxial compression, Computational Materials Science 45 (2009) 321.
(7) S. H. Park, J. S. Kim, J. H. Park, J. S. Lee, Y. K. Choi, O. M. Kwon, Molecular
dynamics study on size-dependent elastic properties of silicon nanocantilevers,
Thin Solid Films 492 (2005) 285.
(8) S. H. Park, H. J. Kim, D. B. Lee, J. S. Lee, Y. K. Choi, O. M. Kwon, Hetero-

68

参考文献

69

geneous crystallization of amorphous silicon expedited by external forcefields: amolecular dynamics study, Superlattices and Microstructures 35 (2004) 205.
(9) Satoshi Izumi and Sidney Yip, Dislocation nucleation from a sharp corner insilicon, Journal of Applied Physics 104, 033513 (2008) .
(10) Ken Gall, M. F. Horstemeyer, Mark Van Schilfgaarde and M. I. Baskes, Atomistic simulations on the tensile debonding of an aluminum-silicon interface, Journal of the Mechanics and Physics of Solids 48 (2000) 2183.
(11) 藤井朋之, 秋庭義明, 単結晶シリコンの破壊挙動に及ぼす欠陥寸法の影響に関する分子動力学解析, 日本機械学会論文集 (A 編) 72-720 (2006),1131.

(12) Yoshitaka Umeno, Takayuki Kitamura, Ab initio simulation on ideal shear
strength of silicon, Materials Science and Engineering B 88 (2002) 79.
(13) 楠幸久, 渋谷陽二, 冨田佳宏, 第一原理計算によるシリコン単結晶の理想せん断
強度, 熱流体系および固体系のミクロシミュレーションに関する合同シンポジウ
ム・分子動力学シンポジウム講演論文集 2000(5) (2000) 53.

(14) 大穂正史，神戸大学修士論文, (2002).
(15) 屋代如月, 北村隆行, 大谷隆一, 野田茂之, 面心立方晶の転位発生の分子動力学シ
ミュレーションおよび前駆局所ひずみ集中の格子不安定条件による特性解析, 日本機械学会論文集，A 編，64-626 (1998), 2583.

(16) 屋代如月, 北村隆行, 大谷隆一, 引張軸に対して直角方向に拘束を受けるニッケル
結晶のへき開破壊に関する原子レベル解析, 日本機械学会論文集, A 編, 63-615(1997), 2380.
(17) 屋代如月, 冨田佳宏, 転位の発生および運動時の局所格子不安定性に関する分子動力学研究, 日本機械学会論文集, A 編, 67-656 (2001), 678.

参考文献

70

(18) Li, J., Van Vliet, K. J., Zhu, T., Yip, S. and Suresh, S., Atomistic Mechanisms
Governing Elastic Limit and Incipient Plasticity in Crystals, Nature, 418 (2002),
307.
(19) Daw, M. S. and Baskes, M. I., Semiempirical Quantum Mechanical Calculationof Hydrogen Embrittlement in Metals, Phys. Rev. Lett., 50–17(1983), 1285.
(20) Daw, M. S. and Baskes, M. I., Embedded–atom method : Derivation and Application to Impurities, Surfaces and Other Defects in Metals, Phys. Rev. B, 29–12(1984), 6443.
(21) Tersoﬀ, J., New empirical approach for the structure and energy of covalentsystems, Phys. Rev. B, 37–12 (1987), 6991.
(22) Tersoﬀ, J., Empirical interatomic potential for silicon with improved elastic
properties, Phys. Rev. B, 38–14 (1988), 9902.
(23) Tersoﬀ, J., Modeling solid–state shemistry : Interatomic potentials for multicomponent systems, Phys. Rev. B, 39–8 (1988), 5566.
(24) Voter, A. F. and Chen, S. P., Accurate Interatomic Potentials for Ni, Al andNi3 Al, Mat. Res. Soc. Symp. Proc., 82 (1987), 175.
(25) Chen, S. P., Voter, A. F., Albers, R. C., Boring, A. M. and Hay, P. J., Investigation of the Eﬀect of Boron on Ni3 Al Grain Boundaries by Atomistic Simulations,J. Mater. Res., 5–5 (1990), 955.
(26) Voter, A. F., Chapter 4 The Embedded Atom Method, John Wiley & Sons,(1994), 77.
(27) Foiles, S. M., Calculation of the Surface Segregation of Ni–Cu Alloys with theUse of the Embedded–Atom Method, Phys. Rev. B, 32–12 (1985), 7685.

参考文献

71

(28) Rose, J. H. Smith, J. R., Guinea, F. and Ferrante, J., Universal features of theequation of state of metals, Phys. Rev. B, 29 (1984), 2963.
(29) Balamane, H., Halicioglu, T. and Tiller, W. A., Phys. Rev. B, 46 (1992), 2250.
(30) 丸山茂夫, 井上知洋, 熱流体系および固体系のミクロシミュレーションに関する
合同シンポジウム・第 5 回分子動力学シンポジウム, (2000), 11.

(31) Smith, R., Nucl. Instr. and Meth., Phys. Res. B, 67, 335 (1992).
(32) Sayed, M., Jeﬀerson, J. M., Walker, A. B. and Cullis, A. G., Molecular dynamics simulations of implantation damage and recovery in semiconductors, Nucl.
Instrum. Meth. Phys. Res. B, 102 (1995), 218.
(33) 上田顯, コンピューターシミュレーション, (1990), 朝倉書店.
(34) Born, M. and Huang, K., Dynamical Theory of Crystal Lattices, (1954), Oxford
UP.
(35) Wang, J., Li, S.Tip, A.R.Phillpot and Wolf, D., Phys. Rev., 52-17B(1995),
12627.
(36) Milstein, F., Hill, R.,Theoretical properties of cubic crystals at arbitrary
pressure-II. Shear moduli, Phys. Rev. B, 3-4B(1971), 1130.
(37) Hill, R., On the Elasticity and Stability of Perfect Crystals at Finite Strain,
Mathematical Proceedings of the Cambridge Philosophical Society, 77 (1975)
225.
(38) Thompson, J. M. T. and Shorrock, P. A., Brifurcaional Instability of An Atomic
Lattice, Journal of Mechanics and Physics of Solids, 23 (1975) 21.
(39) Milstein, F. and Farber, B., Theoretical fcc rightlaw bcc Transition under [100]Tensile Loading, Physical Review Letters, 44-4 (1980), 277.

参考文献

72

(40) Wang, J., S.Tip, S.R.Phillpot and Wolf, D., Phys.Rev.lett., 71-25, (1993), 4182.
(41) Wang, J., Li, S.Tip, A.R.Phillpot and Wolf, D., Phys. Rev., 52-17B(1995),
12627.
(42) Wallace, D. C., Thermodynamics of Crystals, (1972), Wiley, Newyork.
(43) 泉聡志, 東京大学博士論文, (1900).
(44) K. Yashiro, T. Kitamura and R.Ohtani, Atomic Simulation on Deformation
and Fracture of Nano-Single Crystal of Nickel in Tension, JSME International
Journal, Series A, 40-4, (1997), 430.

関連発表講演論文
講演論文

 大石直樹, 屋代如月, 冨田佳宏, 分子動力学法によるシリコンナノカンチレバーの変形・破壊シミュレーション, 関西支部第 84 期定時総会講演会, (2009.3).

 大石直樹, 屋代如月, 分子動力学法による Si 単結晶の理想格子不安定解析, 第 53回日本学術会議材料工学連合講演会, (2009.10).

73

関連発表講演論文
関西-1/174

関連発表講演論文
材料-1/275

関連発表講演論文
材料-2/276

謝辞
本研究を遂行するにあたり，浅学非才の著者が本論文を完成させるまで研究全般に渡り直接御指導頂いたのみならず，日頃の研究生活において，多くの御助言・御協力を頂きました屋代如月准教授に最大限の敬意と感謝の念を表すとともに心より御礼申し上げます．
研究活動に対して幅広い視点からの御助言により，著者の狭い見解を広げて下さった冨田佳宏教授 (現福井工業大学)，阪上隆英教授，長谷部忠司准教授に感謝いたします．また，日々の研究活動が円滑に行えるよう環境づくりに御尽力下さった古宇田
由夫技術職員に心から御礼申し上げます
同じ研究グループの先輩として懇切丁寧に御助言・御支援いただきました西村正臣
さん (現 信州大学)，鈴木雄風さん (現 三菱重工業株式会社) をはじめとする固体力学研究室内諸先輩方に深く感謝致します．また，研究室に配属後約 3 年間，共に研究活動に取り組んだ北村真瑠久君，坂本忠鷹君，富永智之君，平山敬史君，山本隆史君，
山本智君，金谷敬介君に感謝いたします．
最後に，大学生活において積極的に支援して頂いた家族に心から感謝の意を表します．ありがとうございました．

2010 年 2 月
大石77

直樹