修 士 論 文(題 目)

単結晶 Mg ナノワイヤの座屈シミュレーション：
原子弾性剛性係数の固有値・固有ベクトルによる
座屈前駆挙動評価
令和 4 年度

岐阜大学大学院
自然科学技術研究科
物質・ものづくり工学専攻

氏名 太田恵里目次
1

緒言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

解析手法の基礎 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.1

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

4

2.2

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

5

2.3

EAM ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52.4

弾性剛性係数と格子不安定性. . . . . . . . . . . . . . . . . . . . . .

62.5

原子応力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72.6

原子弾性係数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82.7

原子弾性剛性係数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92.8

固有値と固有ベクトル，変形モード . . . . . . . . . . . . . . . . . . .

93

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

114[0001] 圧縮 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

4.1

応力ひずみ曲線，最小固有値および不安定原子数の変化 . . . . . . . . 13

4.2

領域毎の AES の固有値 . . . . . . . . . . . . . . . . . . . . . . . . . . 17

4.3角部原子の前駆挙動

4.4

角部原子の不安定モード . . . . . . . . . . . . . . . . . . . . . . . . . 21

4.5

局所変形発生時の不安定モード . . . . . . . . . . . . . . . . . . . . . 25

. . . . . . . . . . . . . . . . . . . . . . . . . . . 19

[1̄21̄0] 圧縮 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

27

5.1

応力ひずみ曲線，最小固有値および不安定原子数の変化 . . . . . . . . 27

5.2

領域毎の AES の固有値 . . . . . . . . . . . . . . . . . . . . . . . . . . 30

5.3

角部原子の前駆挙動

5.4

角部原子の不安定モード . . . . . . . . . . . . . . . . . . . . . . . . . 33

. . . . . . . . . . . . . . . . . . . . . . . . . . . 31i

[1̄010] 圧縮 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

36

6.1

応力ひずみ曲線，最小固有値および不安定原子数の変化 . . . . . . . . 36

6.2

角部原子の前駆挙動

6.3

角部原子の不安定モード . . . . . . . . . . . . . . . . . . . . . . . . . 40

6.4

領域毎の AES の固有値 . . . . . . . . . . . . . . . . . . . . . . . . . . 42

6.5

領域毎の固有ベクトル成分の変化 . . . . . . . . . . . . . . . . . . . . 43

. . . . . . . . . . . . . . . . . . . . . . . . . . . 39

結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

参考文献 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

謝辞 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

7

ii

1

緒言
一般的なスケールでの材料の変形シミュレーションには連続体モデルである有限要素法などが用いられるが，解析対象がナノスケールになると連続体近似が適用できなくなる．そのため，微視的な材料挙動を解析する最も有力な手法の 1 つとして分子動力学 (molecular dynamics，略して MD) 法が広く用いられている．分子動力学法は各原子についてニュートンの運動方程式 F = ma を立て，数値積分すること
により時々刻々の原子の位置と速度のデータを求める作業を計算機を用いて繰り返
し，物質全体としての性質を導き出す方法である．ここで問題となるのが原子に働く力の評価の精度である．量子力学に基づいて電子状態から力を計算する第一原理
分子動力学法は，計算量が膨大になるためたかだか数千原子に限られる．このためナノスケールの変形，破壊のシミュレーションには原子間に作用する力を近似的に評価する原子間ポテンシャルが通常用いられる．これらのポテンシャルは格子長さや弾性係数，表面エネルギーや昇華エネルギーなど，できるだけ多くの物性を再現するようにフィッティングされてはいるものの，フィッティング範囲から大きくずれた変形，破壊現象においては現実の再現性は保証できず，「似ているが異なる」モデ
ル材料のものとなる．従って，分子動力学法は例えば新しい転位構造の形成過程など「定性的な」局所観察には威力を発揮するが，その発生条件に対する局所エネル
ギーや局所応力などの「定量的な」考察は，現実の材料との違いが必ず存在するた
めあまりなされない．
我々のグループでは，原子シミュレーションにおける定量的な評価の基準として
原子弾性剛性係数 (Atomic Elastic Stiﬀness，以下 AES と称する)Bijα =∆σiα /∆εj に着目した研究を行ってきた (1)–(5) ．ここで ∆σiα は原子応力変化，∆εj はひずみ変化
で，i, j は i, j = 1 ∼ 6 = xx, yy, zz, yz, zx, xy の Voigt 指標 (6), (7) である．6×6 マトリ
クス Bijα の固有方程式 Bijα ∆εj =η α ∆εi が負の固有値を持てば，∆σiα =Bijα ∆εj =η α ∆εiより変形方向と応力方向が逆になる不安定状態と考えられる．対応する固有ベクト

1

ル {∆ε1 , ∆ε2 , · · · , ∆ε6 }={∆εxx , ∆εyy , · · · , ∆γxy } からさらに 3×3 のひずみテンソル[∆εij ] を構築し，その主軸から変形モードを議論してきた．この AES の観点から座屈挙動を議論することを目的として前報 (8) では，hcp Mg 単結晶ナノワイヤの圧縮分
子動力学シミュレーションを行い，S 字状に湾曲した内側部分からの塑性変形（hcp⇒ fcc の相変態）の発生，および，その時の η α(1) <0 の原子の固有ベクトルによる変形モードなどを明らかにしている．
分子動力学法による座屈シミュレーションとして，Wang らは，内圧または外圧を受けた個々の多層カーボンナノチューブ (Multi wall carbon nanotube: MWCNT)の軸方向圧縮座屈を解析した (9) ．Han らは欠陥のある単層カーボンナノチューブ（Single walled carbon nanotube: SWCNT）と二層カーボンナノチューブ (Double

walled carbon nanotube: DWCNT) の座屈を欠陥のないチューブと比較している (10) ．
Ajori らは，ねじり荷重下のホウ素窒化物ナノチューブの機械的特性と座屈挙動を評価している (11) ．Goel らは欠陥のある単層カーボンナノチューブの軸圧縮下での座
屈破壊を観察している (12) ．Wang らは単軸圧縮下の GaN ナノワイヤの座屈特性を調べ，線材の長さが増加するにつれて臨界応力が減少することを報告している (13) ．

Jing らはシリコンナノワイヤの座屈動作に対するシミュレーション温度，ひずみ速度，およびワイヤ長の影響を調査し，温度の上昇とひずみ速度の低下に伴い，臨界荷重が明らかに減少することを報告している．さらに，座屈は動的プロセスの結果として発生するため，材料変形のメカニズムは温度とひずみ速度の両方に影響を受けることを報告している (14) ．
西村らは我々のグループと同様に AES による検討を MWCNT の座屈について行っており，座屈が生じるのは全ての原子の Bijα の最小固有値 η α(1)（西村の論文では λα(1) ）が負になった後，第二固有値 η α(2) も負となる原子が出現したときであること，η α(1)の固有ベクトルは円柱座標で {0, 0, 0, 0, ∆γzr , 0} の面内せん断に対するものであるのに対し，座屈を生じる η α(2) の固有ベクトルは {0, 0, 0, 0, 0, ∆γrθ } の面外せん断に対2するものであることなどを報告している (15), (16) ．さらに，西村らは MWCNT の変形時に生じる局部座屈について軸方向圧縮シミュレーションおよび曲げ変形シミュ
レーションを実施し，AES の観点から検討している．その結果，η α(2) <0 となる原子が多くの層に発生した領域から局部座屈を生じること，曲げ変形解析で MWCNT
全体の曲げモーメントの変化が線形から非線形になるとき，η α(2) が負となる原子が急増することを明らかにしている (17) ．
前報では局所変形開始に焦点をあてて議論を行ったが，局所変形開始以前に系はS 字状に湾曲していた．そこで本研究では前報と同様に hcp Mg 単結晶細線の圧縮シ
ミュレーションを結晶方位を変えて行い，細線が湾曲し始める点について η α(1) < 0の原子，すべての原子の固有値・固有ベクトルから検討する．

3

2

解析手法の基礎2.1

分子動力学法

分子動力学法は，系を構成する各原子についてニュートンの運動方程式mαd2 r α
= Fα
2
dt

(2.1)

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

Fα = −

∂Etot
∂r α

(2.2)

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

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

(2.3)
(2.4)

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

dr α= v α (t)
dt

(2.5)

であり，式 (2.1) と式 (2.5) を式 (2.3) と式 (2.4) に代入すると，()(∆t)2 F α (t)
3
+
O
(∆t)
2
mα
2
(
)
(∆t) F α (t)3r α (t − ∆t) = r α (t) − ∆tv α (t) +
+
O
(∆t)
2
mα

r α (t + ∆t) = r α (t) + ∆tv α (t) +

(2.6)
(2.7)

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

r α (t + ∆t) + r α (t − ∆t) = 2r α (t) + (∆t)2()
F α (t)
4
+
O
(∆t)
mα

(

r α (t + ∆t) − r α (t − ∆t) = 2∆tv α (t) + O (∆t)3
4

)

(2.8)
(2.9)

が得られる．∆t3 以上の高次項は無視できるとすると，時刻 t + ∆t での位置ベクトルと t での速度は
2 F

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

v α (t) =αα

1{r α (t + ∆t) − r α (t − ∆t)}
2∆t

α

(t)

mα

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

2.2

原子間ポテンシャル

式 (2.2) で示したように，原子 α に作用する力 F α は系のエネルギー Etot をポテン
シャルとして決定される．したがって，系のポテンシャルエネルギー Etot をいかに精度よく評価するかが重要となる．量子力学に基づき，電子や原子核のハミルトニ
アンから系のポテンシャルエネルギーを精密に求めて原子の運動を追跡する第一原
理分子動力学 (20) も試みられているが，計算量が極めて膨大になるため，ごく少数の原子しか扱うことができず，変形・破壊のような多数の原子の動的挙動への直接的な適用は困難である．そこで，原子間相互作用を簡略評価する原子間ポテンシャ
ルが通常用いられる．

2.3

EAM ポテンシャル

EAM(Embedded Atom Method) は金属中の多体効果を良好に再現することから広く用いられている．密度汎関数理論に基づき，まず金属材料における系のポテン
シャルエネルギー Etot は原子を価電子雲中に埋め込むエネルギーと原子間の 2 体間相互作用の和で与えられるとする．さらに，埋め込みエネルギーは埋め込む位置の5電子密度にのみ依存すると仮定することによって，系全体のエネルギーは次式のように表される．

Etot =∑F α (ρα ) +α1 ∑ ∑ αβ αβϕ (r )
2 α β(̸=α)(2.12)

ここで，ρα は原子 α の位置における多体効果を考慮する密度を表し，F α (ρα ) は密度 ρα の位置に原子を埋め込むエネルギー，ϕαβ (rαβ ) は距離 rαβ 離れた原子 α と β のクーロン相互作用である．密度 ρα は周囲（neighbor）の原子 β からの寄与 f β (rαβ )の重ね合わせで与えられると仮定しρα =neighbor∑f β (rαβ )(2.13)

β(̸=α)

で評価する．
本研究では，Zhou らが多数の元素に対してパラメータを提案しているポテンシャ
ルを用いた (21)–(23) ．

2.4

弾性剛性係数と格子不安定性

応力 σij および弾性係数 Cijkl は，等温過程では

σij

1
=
V

(

∂F
∂ϵij

)

1
, Cijkl =
V

(

∂ 2F
∂ϵij ∂ϵkl)(2.14)と定義される (24) ．ここで，F は Helmholtz の自由エネルギー (断熱過程では内部エネルギー U )，V は結晶の体積，ϵij は平衡状態 (無負荷とは限らない) からの仮想的な微小ひずみである．一方，無負荷平衡状態を基準とするひずみ εij と応力 σij の関
係は，2 つの平衡状態間の変形を考えて導出される次の弾性剛性係数によって表される (24) ．

Bijkl ≡

(

∆σij
∆εkl

)

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

ここで，δij はクロネッカーのデルタである．Wang, Yip らは，Voigt 対称性をもたsym≡ (Bijkl + Blkji )/2 の正値性によって結晶の安定性を評価するこ
せたテンソル Bijklとを提案している (6), (7) ．

2.5

原子応力

簡単のため，結晶の内部エネルギー U が系全体のポテンシャルエネルギー Etot に等しいとする．このとき，応力は平衡状態からの微小ひずみ ϵ に対するポテンシャ
ルエネルギーの単位体積当たりの変化として与えられる (24) ．

σij =1 ∂EtotV ∂ϵij(2.16)

ここで，V は平衡状態における系の体積であり，下付添字の i, j はテンソルのデカ
ルト座標成分を表す．(2.16) 式の微分を求めるため，平衡状態からの仮想的な均一変形を考える．結晶内の原子 α の位置ベクトルは仮想変形のヤコビ行列 J によってr α = J r̄ α

(2.17)と変化する．ここで，「 ¯ 」は仮想ひずみによる変形前の値を示す．これより，原子α と 原子 β の間の距離 rαβ には(rαβ )2 = r̄iαβ Gij r̄jαβ(2.18)なる関係が成立する．ただし，Gij = Jki Jkj である．仮想変形の Lagrange ひずみテンソル ϵij は

ϵij =

]
1[
Gij − δij
2

(2.19)であり，その微小量dϵij =

1dGij
2(2.20)

7

と式 (2.18) の関係から次の関係が得られる．

r̄iαβ r̄jαβ
=
rαβ

∂rαβ
∂ϵij

(2.21)

これより EAM ポテンシャルにおける応力は次式で評価される．
(

N ∑
N
1 1∑
1 ∂Etot∂rαβ ∂Etot=
V ∂ϵijV 2 α β(̸=α) ∂ϵij ∂rαβ

σij =

)

) αβ αβN ∑N (r r
1 ′ αβ αβ1 ∑′α α′ β αβ=ϕ (r ) + F (ρ )f (r ) i αβjV α β(̸=α) 2r

ここで，原子体積を Ω = V /N とすれば，各原子のエネルギー寄与

Etot =∑F α (ρα ) +α1 ∑ ∑ αβ αβϕ (r )
2 α β(̸=α)の微分から次のように原子応力を定義できる．

σijα =

2.6) αβ αβN (r r
1 ∑1 ′ αβ αβ′α α′ β αβϕ (r ) + F (ρ )f (r ) i αβjΩ β(̸=α) 2r

原子弾性係数弾性係数も応力と同様に U ≈ Etot の場合にはCijkl =

1 ∂ 2 Etot
V ∂ϵij ∂ϵkl

(2.22)であるので，平衡状態からの仮想均一変形を考えると EAM ポテンシャルにおける弾性係数は以下のようになる．

Cijkl =(N ∑} αβ αβ )
N ∑
N
N {
ri rj
1 ∑
∂rαβ ∂ ∑1 ′ αβ αβ′α α′ β αβ
ϕ
(r
)
+
F
(ρ
)f
(r
)
2V α β(̸=α) ∂ϵkl ∂rαβ α β(̸=α) 2rαβ[} αβ αβ αβ αβ
N ∑
N {
1 1∑
ϕ′ αβ (rαβ ) ri rj rk rl=
ϕ′′ αβ (rαβ ) −V 2 α β(̸=α)
rαβ(rαβ )2

+

N ∑
N
∑

′α

{
α

F (ρ ) f′′ β

α β(̸=α)

+

N
∑
α

F

′′ α

{ ∑

} αβ αβ αβ αβ
f ′ β (rαβ ) ri rj rk rl(r ) −rαβ(rαβ )2αβ]

riαβ rjαβ }{ ∑ ′ β αβ rkαβ rlαβ }(ρ )f (r ) αβf (r ) αβrrβ(̸=α)
β(̸=α)
α

′β

αβ

8

(2.23)

応力と同様に，各原子位置における原子弾性係数を以下のように定義する．
[

α
Cijkl

} αβ αβ αβ αβ
N {ϕ′ αβ (rαβ ) ri rj rk rl
1 1 ∑ϕ′′ αβ (rαβ ) −=Ω 2 β(̸=α)
rαβ(rαβ )2

+

N
∑

′α

{
α

F (ρ ) f′′ ββ(̸=α)} αβ αβ αβ αβ
f ′ β (rαβ ) ri rj rk rl(r ) −rαβ(rαβ )2αβ{ ∑αβ αβ }{αβ αβ }∑′′ α α′ β αβ ri rj′ β αβ rk rl+ F (ρ )f (r ) αβf (r ) αβrrβ(̸=α)
β(̸=α)](2.24)

原子弾性剛性係数

2.7以上で定義した原子応力，弾性係数から，Voigt 対称性をもたせた原子弾性剛性係数を以下で求める．
α
α
α
α
α
Bijkl
= Cijkl+ (σilα δjk + σjl
δik + σik
δjl + σjk
δil − 2σijα δkl )/2(2.25)αBijkl
は 4 階のテンソルであるが，Voigt 対称性を持たせた場合，弾性係数テンソルと同じく独立な成分は 21 個となり，xx，yy，zz，yz，zx，xy を 1 ∼ 6 とする Voigt
表記を用いれば 6 × 6 のマトリックスとして Bijα と表すことができる．その場合の式(2.25) は以下のようにマトリクスで表すことができる．

[

]

[

Bijα = Cijα

]



1
+ 
2



2σ1α

−σ1α − σ2α2σ2α−σ1α − σ3α−σ2α − σ3α
2σ3α

sym.

−σ4α
σ4α
σ4α
α
σ2 + σ3α

σ5α
−σ5α
σ5α
σ6α
α
σ3 + σ1α


σ6α

σ6α

α

−σ6
(2.26)
α

σ5

α

σ4
σ1α + σ2α固有値と固有ベクトル，変形モード

2.8

固有方程式
[




∆ε1 





]

Bijα 

..

. 




 ∆ε 
6




∆ε1 







= ηα 

..

(2.27).  ∆ε 
6

9

を解くことは Bijα を対角化することと等価である．

[

Bijα

]




⇒




η α(1)000

η α(2)
η

α(3)

0
0
0

0
0
0
0

η α(4)η α(5)

0
0
0
0
0










(2.28)

η α(6)

sym.

せん断成分は互いに独立なので，6×6 の Bijα マトリクスは通常は次のような形となる．


α
B11
α
 B21
 α
[
]

B33
Bijα = 
 0

 0
0

α
B12
α
B22
α
B32
0
0
0

α
B13
α
B23
α
B33
0
0
0

0
0
0
α
B44
0
0

0
0
0
0
α
B55
0


0
0 

0 

0 

0 
α
B66

(2.29)αα
したがって B44
∼ B66
は固有方程式 B αij ∆εj =η α ∆εi の解で，それぞれの固有ベクト

ルは {0, 0, 0, ∆γyz , 0, 0}，{0, 0, 0, 0, ∆γzx , 0}，{0, 0, 0, 0, 0, ∆γxy } である．垂直ひずみに対する 3×3 の部分マトリクスの固有値は，無ひずみの等方弾性体では特性方程式



det 


C11 − ηsym.

C12C11 − η



C12

C12 
=0
C11 − η(2.30)の解から η=C11 −C12 (重解)，C11 +2C12 となる (C11 , C12 は弾性係数)．この C11 +2C12に相当する固有値は他の固有値より大きく (第 6 固有値)，負になることはない．等αα方でない場合，C11 − C12 に相当する固有値は重解とならず，B44
∼ B66
との大小関

係で第 1∼5 固有値のいずれかとなる．(2.27) 式の固有ベクトルの ∆ε4 ∼ ∆ε6 は工学
せん断ひずみ ∆γyz ，∆γzx ，∆γxy なので，1/2 して

∆εxx ∆εxy ∆εxz
]
[



∆εij = ∆εyy ∆εyz 

sym.
∆εzz

(2.31)の 3×3 テンソルとし，これの主軸，すなわちこのテンソルの固有ベクトル x1 , x2 , x3を求めることができる (固有値 ∆ε1 < ∆ε2 < ∆ε3 の順)．本論文では最大せん断方向の x1 + x3 を変形モードとして定義する．

10

3

シミュレーション条件
図 3.1 に示すような xy 断面が正方形の細い角柱状の Mg 単結晶を解析対象とする．

また，表 3.1 に解析を行った各モデルの結晶方位，外形，原子数を示す．結晶方位はそれぞれ x,y,z を hcp の [1̄21̄0],[1̄010],[0001] とする座標系，[0001],[101̄0],[1̄21̄0] とする座標系，[0001],[1̄21̄0],[1̄010] とする座標系の 3 通りである．できるだけ断面寸法，z 方向長さ，原子数が同じになるように図 3.1(b) に示す単位格子を並べた．原子間相互作用の評価には Zhou(21) らの EAM ポテンシャルを用いた．また，すべてのシミュ
レーションにおいて全方向周期境界下で真空中で孤立した細線を模擬するため，シミュレーション空間の x,y 方向は細線断面より十分大きくしている．熱による影響を排除するため温度は 0.1[K] とし，速度スケーリング法により制御した． 20,000[step]の緩和計算を行った後，z 方向に毎ステップ ∆εzz =1.0 × 10−7 の割合で原子座標・セル辺長をスケーリングして圧縮シミュレーションを行った．

Table 3.1

Simulation conditions.

11Fig. 3.1

Simulation model [0001].

12

4[0001] 圧縮

4.1

応力ひずみ曲線，最小固有値および不安定原子数の変化

図 4.1 に η α(1) の最小値の変化 (緑線)，図 4.2 に η α(1) <0 の原子数の変化 (青線) を，
応力―ひずみ曲線 (黒線) と圧縮された角柱の湾曲度合い（オレンジ線）の上に重ねて示した．湾曲度合いは角柱の軸方向の各原子面（角柱断面）の重心位置 (xi , yi ) を，
ひずみ 0 における位置 (xi0 , yi0 ) からの変位 ∆ri =√(xi − xi0 )2 + (yi − yi0 )2 として評価し，各ひずみにおいて ∆r の最大値（最も軸からずれた原子面の値）をプロットしている．赤い縦線は塑性変形が確認されたひずみを示し，左軸目盛が応力，右軸目
盛内側が角柱断面中心の最大ずれである．右軸目盛り外側に最小固有値，η α(1) <0 の原子の数を示している．ひずみ，応力は圧縮を正としている．応力は圧縮ひずみとともに単調に増加し，局所変形=塑性変形を生じ始めると応力ピーク σzz =2.32[GPa]を示し応力は急減する（図中 (f)）．オレンジ線で示す湾曲度合いを見ると，ひずみ0.03 付近から急激に湾曲し始めていることが分かる．このとき，応力勾配も減少しα(1)非線形性を示す．図 4.1 の ηmin はすべての原子の中の最小値を示しているので同一α(1)原子の η α(1) を追跡したものではないが，この ηmin は圧縮前から負値を示しており，
α(1)後の図 4.3 で示すように角部の原子のものである．ηmin は圧縮によってひずみ 0.02までさらに低下し，(b) 点付近からやや回復してほぼ一定になり，局所変形を生じる(e) 点で急減する．図 4.2 の η α(1) <0 の原子数の変化について，図中 (a)∼(f) で示した
ひずみでのスナップショットを η α(1) の正負で色分けして図 4.3 に示した．η α(1) < 0
原子を赤，η α(1) > 0 原子を緑で着色している．上図 (i) は xz 側面を，下図 (ii) は xy
断面を上から見た図である．図 4.3(ii)(a) で確認できるように，圧縮前から角柱の四隅の 1 原子列が η α(1) < 0 となっており，この総数が 800 個で図 4.2 で (d) 点までその数が変わらない．図 4.3(i)(d) で見た目にも S 字状になり始めた角柱の圧縮側角部に η α(1) < 0 の原子が発生する（図中矢印）．その後，図 4.3(i)(f) の矢印で示す箇所

13

に局所変形を生じると η α(1) < 0 原子は急増した．なお，(ii) の断面で (a)∼(f) いず
れにも斜めのずれが見えるが，周期境界下で表示されている原子の段であり塑性変形による結晶のずれではない．

14

Fig. 4.1

Stress-strain curve，deviation from axis center and change in the minimum of the 1st eigenvalue η α(1) under [0001] compression.

Fig. 4.2

Stress-strain curve，deviation from axis center and the number of η α(1) <0 atoms under [0001] compression.

15

(a)

(b)

(c)

(d)

(e)

(f)

(i) Side view

(a)

(b)

(c)

(d)

(e)

(f)

(ii) Top view (Cross section under periodic boundary in z-axis.)
Fig. 4.3

Snapshots of Mg nanowire at Points (a)∼(f) in Figs. 4.1 and 4.2. Atoms
are colored by the sign of η α(1) , green ; η α(1) > 0, red ; η α(1) < 0.

16

4.2

領域毎の AES の固有値

図 4.4 は図中右の断面で色分けした場所ごとに η α(1) を平均して圧縮時の変化を示したグラフである．図には先述の角柱断面中心の最大ずれをオレンジ線で示している．圧縮前の値は角部（corner）の η α(1) が最も低く，続いて x 方向表面（(1̄21̄0) 面），

y 方向表面（(1̄010) 面），内側のバルク部の順に高くなる．なお角部は頂点の 1 原子
列は η α(1) < 0 であるが，1 格子分（4 列）の平均をとると正値となっている．引張を与えた場合は固有値が減少するが (25) ，[0001] 圧縮では corner 以外は η α(1) の値がα(1)増加する．先の図 4.1 では ηmin は低下していたので，角部も頂点を除く 3 原子列は上昇している．すなわち，角柱の角部の 1 原子列のみ図 4.1 で示したように固有値は低下する．圧縮初期は特に角部に見られるように振動しているが，ひずみ 0.02 以
降は収まっている．また，角部以外の η α(1) は振動が収まるひずみ 0.02 付近が頭打ちになる．まずバルク部で η α(1) の上昇が頭打ちになり，少し遅れて (1̄21̄0)，(1̄010)
表面部が続く．この後にずれが急増しはじめる．赤い破線以降は内部の構造変化を生じて η α(1) に揺らぎを生じる．
図 4.5 に完全結晶の静力学解析で，横方向応力０条件で [0001] 圧縮したときの第 1∼
第 5 固有値の変化を示す．x[1̄21̄0], y[1̄010],z[0001] の座標系において，無ひずみでの√2
A1 − A22 +8C13
固有値は η α(1),(2) = C44 , η α(3) = C66 , η α(4) = C11 − C12 , η α(5) =, (A1 =2
√ 2 2
A1 + A2 +8C13
C11 + C12 + C33 , A2 = C11 + C12 − C33 ) となる．η α(6) は η α(5) の共役解 ()
2で常に正の大きな値を示すため省略している．静力学解析ではひずみ 0.03 手前で C44（zx，yz 面のせん断）から C66（xy 面のせん断）に固有値の大小関係がスイッチしており，図 4.4 で角部以外の固有値が頭打ちになるのはこれに対応していると考える．

17

Eigenvalue, η (1)ave, GPa

er

Inn

15

20



Conner

10

5

Maximum deviation of cross-sectional centerfrom initial position, r, Å

30

20

10

Corner

0

0.02

0.04

0
0.06

Inner

Strain, εzz(a) Average of η α(1) for each zone in (b) (b) Zone in cross section for η α(1) evaluaand deviation from axis center.
tion.

Fig. 4.4

Change in η α(1) averaged by conner, surface and inner under [0001] compression.

Fig. 4.5

Change in the eigenvalues of Mg bulk perfect lattice under [0001] compression (static analysis).

18角部原子の前駆挙動

4.3

図 4.6 に η α(1) < 0 原子の固有ベクトルの成分の絶対値 |∆εxx |, |∆εyy |, · · · , |∆γxy | の平均の変化を示す．図には湾曲度合いを示すグラフ（右目盛り）と塑性変形を生じたひずみ，先の図 4.1∼4.3 の (a)∼(f) 点をあわせて示している．固有ベクトルは大きさが√∆εxx 2 + · · · + ∆γxy 2 = 1 となるように規格化されているので値の大きさに意

味はなく，また正負逆方向をとりうるため絶対値での平均とした．湾曲度合いが増えはじめる (b) 点で成分の大小関係の急激な入れ替わり＝変形モード変化を生じていることが分かる．ひずみ 0.02 までは |∆γxy | が最も高い値を示し，|∆εyy | は 0，他
は 0.2∼0.4 の値をとっている．ひずみ 0.02 以降は |∆γxy | はさらに増加するが，それ以上に顕著な変化として z 方向への垂直ひずみ |∆εzz | と zx 面でのせん断 |∆γzx | の急
増，yz 面のせん断 |∆γyz | の急減を生じている．|∆εyy | 方向へのモードもひずみ 0.02以降に増加している．図 4.3 の側面図で S 字状になりはじめるのは (c) → (d) 間であるが，ここで yz 面でのせん断 |∆γyz | が著しく増加し，|∆γzx | と |∆εzz | が急減する．
図 4.3(ii)(d)∼(f) ではらせんのように変形しているが，zx のせん断から yz のせん断にシフトしたものと考えられる．

19

Eigenvector,|γyz||γzx||γxy|(d)(e)(f)

0.6

30

0.4

20

0.2

10

0

0.02

0.04

Strain, εzz
Fig. 4.6

0
0.06

Maximum deviation of cross-sectional centerfrom initial position, r, Å(a)|εxx||εyy||εzz|(c)(b)

Change in the components of the 1st eigenvector of η α(1) < 0 atoms anddeviation from the axis center under [0001] compression. Points (a)∼(f)correspond to those in Figs. 4.1∼4.3.

20

4.4

角部原子の不安定モード

図 4.3 で示したスナップショットについて，角部の η α(1) < 0 原子の固有ベクトル{∆ε1 , ∆ε2 , · · · , ∆ε6 }={∆εxx , ∆εyy , · · · , ∆γxy } から求めた最大せん断方向 (8) を，固有値の大きさ |η α(1) | でスケーリングした矢印で可視化した図 (不安定モード) を図 4.7に示す．グラフ中の矢印が図 4.3 の赤原子の不安定モード，緑の丸は赤原子の位置を表す．図 (ii) xy 断面の黒枠は見やすくするためにつけたものであり，周期境界ではないことに注意されたい．左側表面角部の原子は右側表面のそれより矢印が著しく大きいが，これは後で説明するように単位格子の位置関係による．圧縮前の図 (a)を見ると，左側表面の角部の原子はほぼ 2 方向に向いている．後述するように (a) →(b) では左側の η α(1) < 0 の原子列の一部の原子が内側にずれ，方向も 4 方向ぐらいになる．(b) → (c) は前節でモード変化を生じた後であり，角部の不安定モードは xy
断面で見ると放射状に様々な方向を向いている．その後角柱は S 字状に湾曲し始め，
図 (i)(e),(f) では新たに η α(1) < 0 となった部分が見られるが，特に大きな矢印は認められない．なお，xy 断面で楕円状に見えるのは図 4.3(ii) のようにらせん状になった角柱の赤い原子を全て投影しているためである．

21

(a)

(b)

(c)

(d)

(e)

(f)

(i) Side view

(a)

(b)

(c)

(d)

(e)

(f)

(ii) Top view (Cross section under periodic boundary in z-axis.)
Fig. 4.7

Unstable mode under [0001] compression. Maximum shear direction (unstable mode) calculated from the eigenvector of {∆εxx , ∆εyy , · · · , ∆γxy }for η α(1) < 0 atoms in Fig.4.3(red atoms). The position of red atoms
are indicated by green dot. The shear directions are projected on xz andxy planes by purple arrows. The length of the arrow are scaled by themagnitude of |η α(1) |.
22

図 4.8 に図 4.3(a),(b) の左側表面の一部を拡大して示す．(b) の拡大図の破線で示した丸は，ずれる前の η α(1) < 0 原子の位置である．(a) では η α(1) < 0 原子は一列に並んでいるが，(b) では図中の四角で囲った部分のように η α(1) < 0 の原子列の一部の原子が原子半個分内側にずれている．図 4.9 は図 4.3(b) における不安定モードを拡大して示したものである．角部の η α(1) < 0 原子の 1∼4 の原子列は左右の模式図に示した単位格子の 1∼4 の位置に対応する．単位格子上には η α(1) < 0 原子の不安定モードの方向を模式的に示した．左側表面角部の原子は右側表面のそれに比べて「飛び出た」状態になっているため構造的に不安定である．このため |η α(1) | が著しく
大きく，大きな矢印の不安定モードとなる．また 3，4 の不安定モードの方向は空き
サイトへの方向を示しているが，1，2 の原子列は先の図 4.8 の原子移動もあり必ず
しも隣接サイト方向を向いていない．

23

(a)
Fig. 4.8(b)Migration of corner atoms during Fig. 4.3(a) → (b).

Fig. 4.9

Details of unstable mode in Fig. 4.3(b).

24

4.5

局所変形発生時の不安定モード

図 4.10 に応力急減時の局所変形発生部（図 4.3(i)(e) の四角で囲った部分）を拡大
したスナップショット（原子は η α(1) の正負で着色）に不安定モードを重ねて示す．

(a) は xz 側面，(b) は A-A′ 断面を上からみたものである．不安定モードは (a) は xz 表面のみ，(b) では (a) の A-A′ より下部分について透過して示している．(1̄010) 表面の１原子列内側に η α(1) < 0 の原子が A-A′ の位置に発生して縦に広がる（図 (iii)）．
圧縮とともにさらに内側へ η α(1) < 0 原子列が発生し，そこを起点として局所変形を生じた（図 (v)，A-A′ 断面左下角部）．この領域では結晶構造が hcp から fcc に変化しており，境界部には η α(1) < 0 原子が見られる．次に不安定モードに着目すると，(ii) の不安定モードは後に変形を生じる湾曲凹側で A-A′ 断面を境にモードの向きが異なっている．(iii) で 1 原子列内側に発生した不安定モードの大きさはこの時点では小さい．(iv) でさらに内側に η α(1) < 0 の原子列が生じたときにはこの部分のモードが大きくなり方向も様々である．結晶にずれを生じるときは原子が相対的に様々な方向へ移動するためこのように複雑なモードとなる．(v) は相変化が進行している
状態で，原子移動方向に相当するモードが変形の先端境界部分の η α(1) < 0 原子に見
られる．

25

A

A

(i) εzz =0.0 (ii) εzz =0.040 (iii) εzz =0.041 (iv) εzz =0.042 (v) εzz =0.043(a) Side view(i) εzz =0.0 (ii) εzz =0.040 (iii) εzz =0.041 (iv) εzz =0.042 (v) εzz =0.043(b) A-A′ cross section
Fig. 4.10

Unstable modes of η α(1) <0 atoms in hcp → fcc transformation under[0001] compression.

26

5[1̄21̄0] 圧縮

5.1

応力ひずみ曲線，最小固有値および不安定原子数の変化

図 5.1，図 5.2 に η α(1) の最小値の変化 (緑線) と η α(1) <0 の原子数の変化 (青線) を，
応力―ひずみ曲線 (黒線)，湾曲度合いのグラフに重ねて示す．本圧縮方向では応力
ひずみは正弦波状を示しているが，塑性変形を生じるのは (f) 点である．ひずみ制
御下の圧縮のため，応力ピーク後も局所変形を生じずに応力が減少しているが，(d)α(1)

点で圧縮方向の変形抵抗はすでに喪失している．前章同様，ηmin は圧縮前から負値であり，図 5.1 で (c) 点まで η α(1) < 0 の原子数は 648 個で一定ある．図中 (a)∼(f) で示したひずみでのスナップショットを η α(1) の正負で色分けして図 5.3 に示した．図5.3(ii)(a) に拡大して示したように，圧縮前から角柱の角部の 1 原子列とその対角にある角部が η α(1) < 0 であり総数は 648 個である．図 5.1，5.2 でずれが生じはじめるのは (b) 点であるが，このときは η α(1) < 0 原子の数に変化はない．応力ピークのわ
ずかに前の (c) 点で η α(1) < 0 原子がそれまで正値であった角部に出現し（図 5.3(c)矢印），上下に広がっている（図 5.3(d) 矢印）．ただし，図 5.3(e) に丸印で囲ったよ
うにこの部分の η α(1) は回復し，四角で囲った部分の角部から内部に η α(1) < 0 の原子が発生し (f) の矢印部分に局所変形を生じた．前章では S 字状に湾曲した凹部（圧縮）から局所変形を生じたが，本章の圧縮では凸部（引張）から発生している．ただし，ひずみ制御下での (d)∼(f) の変化は応力制御下で再現される保障はなく，系の不安定挙動の議論は (d) 点までに限られる．

27

Stress, σ , GPa

(d)
2

100

(c)

z
z



(b)
1

50

(e)
(a)
0

0.05

0
0.15

0.1

Strain, εFig. 5.1

-0.1

-0.2

-0.3

-0.4

Stress-strain curve，deviation from axis center and change in the minimum of eigenvalue under [1̄21̄0] compression.

Stress, zz
Maximum deviation of cross-sectional centerfrom initial position, r
Number of unstable atoms, Ninst
150

Stress, σ , GPa

(d)
2

100

(c)

z
z

2500

2000

1500



1

(b)

50

(e)
(a)
0

0.05

0.1

Strain, ε0
0.15

Number of unstable atoms, Ninst(f)

Maximum deviation of cross-sectional centerfrom initial position, r, Å

3Fig. 5.20Minimum eigenvalue, η (1)min, GPa(f)Maximum deviation of cross-sectional centerfrom initial position, r, Å

3Stress, zz
Maximum deviation of cross-sectional centerfrom initial position, r
Minimum eigenvalue, η (1)min
150

1000

500

0

Stress-strain curve，deviation from axis center and the number of η α(1) <0 atoms under [1̄21̄0] compression.

28

(a)

(b)

(c)

(d)

(e)

(f)

(i) Side view

(a)

(b)

(c)

(d)

(e)

(f)

(ii) Top view (Cross section under periodic boundary in z-axis.)
Fig. 5.3

Snapshots of Mg nanowire at Points (a)∼(f) in Figs. 5.1 and 5.2. Atoms
are colored by the sign of η α(1) , green ; η α(1) > 0, red ; η α(1) < 0.

29

領域毎の AES の固有値5.2

図 5.4 は図 4.4 と同様に，図中右の断面で色分けした場所ごとに η α(1) を平均して圧縮時の変化を示した．圧縮前の x，y 表面（(0001) と (101̄0)）の η α(1) の差が小さ
く，圧縮初期の振動が少ない．[0001] 圧縮と異なり，引張で見られるように (25) 固有
値はひずみによって減少している．湾曲が急増しはじめた後，ひずみ 0.05 くらいで極小値を示して増加に転じている．

Eigenvalue, η (1)ave, GPa

150

15

Inn

100

er



Conner

10

5

Maximum deviation of cross-sectional centerfrom initial position, r, Å

20

50

Corn

0

er

0.05

0.1

0
0.15

Strain, εzz(a) Average of η α(1) for each zone in (b) (b) Zone in cross section for η α(1) evaluaand deviation from axis center.
tion.

Fig. 5.4

Change in η α(1) averaged by conner, surface and inner under [1̄21̄0] compression.

30

5.3

角部原子の前駆挙動

図 5.5 に η α(1) < 0 原子の対応する固有ベクトルの成分の絶対値の平均の変化を示す．先述のように，系の不安定挙動を議論できるのは (d) 点までなのでそれ以降は参考とし薄墨をかけている．やはり，ずれが発生する (b) 点で角部の急激なモード
変化を生じている．[0001] 圧縮では湾曲が増加するまでは角部のモード変化は少な
かったのに対し，[1̄21̄0] 圧縮では |∆γyz |，|∆γzx | が湾曲前から著しく変化している．
始めは |∆γzx | の値が一番大きく，次いで |∆γyz | であったが |∆γyz | が急増，|∆γzx | が急減してすぐに入れ替わる．(b) 点では逆に |∆γyz | が急減，|∆γzx | が急増する．図5.3(ii) でも角柱はらせん状に変形しているが yz せん断から zx せん断へのモード変化
が本圧縮での座屈の前駆挙動と考える．

31

(b)

(e) (f)

1

150

Eigenvector,

0.8
100
0.6



0.4
50
0.2

0

0.05

0.1

Strain, εzz
Fig. 5.5

Maximum deviation of cross-sectional centerfrom initial position, r, Å(a)| γyz|| γzx|| γxy|| εxx|| εyy|| εzz|(c)(d)

00.15

Change in the components of the 1st eigenvector of η α(1) < 0 atoms anddeviation from the axis center under [1̄21̄0] compression. Points (a)∼(f)correspond to those in Figs. 5.1∼5.3.

32

5.4

角部原子の不安定モード

図 5.6 に図 5.3(a)∼(c) のひずみにおけるスナップショットに不安定モードを重ねて示す．図は上から順に図 5.3 の四角で囲った部分の (i)yz 表面，(ii)xz 表面，(iii)は xy 断面を上から見たものとなっている．本圧縮では図 4.8 で示したような角部のη α(1) < 0 の原子列のずれは生じていない．圧縮前は yz 面では斜めのモード，xz 面では 4 つは+x 方向を向いているがそれ以外は斜めのモードとなっている．そしてこれ
ら斜めのモードはそれぞれ yz 面，xz 面に沿っているので，直角方向から見ると点に見える．xy 断面から見たモードは 6 方向を向いており上下で対称となっている．(b)は前節で固有ベクトルに変化を生じた点であり，(a) → (b) で圧縮前は+x 方向を向
いていたモードも斜めになる．xy 断面を見ると依然としてモードは 6 方向を向いているが，その方向は圧縮前より等方的になっている．また，+x 方向の矢印が少し短くなる．(b) → (c) で xz 面のモードが見えなくなり，ほぼ yz 面に平行なモードとなる．(b) では隣り合う原子のモードが同じ方向を向いているとは限らなかったが，(c)では隣り合う原子のモードは同じ方向を向いており，ある原子を境にその方向が変
わる．

33

(a)

(b)
(i) Side view of yz plane

(c)

(a)

(b)
(ii) Side view of xz plane

(c)

(a)
(b)
(c)
(iii) Top view (Cross section under periodic boundary in z-axis.)
Fig. 5.6

Distribution of η α(1) <0 atoms and unstable mode in Fig. 5.3.
34

図 5.7 に図 5.3(b) における不安定モードと単位格子の模式図を示した．角部のη α(1) < 0 の原子列は左右の単位格子の位置に対応する．角部の原子 1，2 は他の角部に比べて飛び出た状態になっている．また，[0001] 圧縮では角部によってモードの大きさに違いが見られたが本圧縮では見られない．これは x,y,z が [0001],[101̄0],[1̄21̄0]の結晶方位において原子 1，2 は等価であるからだと考える．

2
2

1

1
Fig. 5.7

Details of unstable mode in Fig. 5.3(b).

356[1̄010] 圧縮

6.1

応力ひずみ曲線，最小固有値および不安定原子数の変化

図 6.1，図 6.2 に η α(1) の最小値の変化 (緑線) と η α(1) <0 の原子数の変化 (青線) を，
応力―ひずみ曲線 (黒線)，湾曲度合いのグラフに重ねて示す．また，図中 (a)∼(f) で示したひずみでのスナップショットを η α(1) の正負で色分けして図 6.3 に示した．応力は圧縮ひずみとともに単調に増加し，応力ピーク σzz =1.53[GPa] を示し降伏する．
α(1)

前章同様，ηmin は圧縮前から負値である．はじめは緩やかに減少していたが，局所変形を生じる直前に急減し変形後は増減を繰り返す．図 6.3(ii)(a) に見られるように，
圧縮前から角柱の四隅の 1 原子列が η α(1) < 0 であり総数は 748 個である．(b) 点は角部の η α(1) < 0 原子の固有ベクトル成分に顕著な変化が生じた点（後述），(c) 点は湾曲開始点であるが，この時点では η α(1) < 0 原子の数に変化はない．η α(1) <0 の原子は応力ピーク前の (d) 点で図 6.3(i)(d) の矢印で示す湾曲圧縮側の角部で増加し始め，
応力ピークの (e) 点では図 6.3(i)(e) の楕円で囲ったように yz 面に η α(1) <0 の原子が広がる．また，図では反対側であるため見えないが，裏側の湾曲圧縮側の角部でも

η α(1) < 0 の原子が同様に発生している．局所変形を生じるわずかに前で η α(1) < 0 原子が急増し，図 6.3(i)(f) の矢印で示す箇所（η α(1) <0 の原子が初めに増加した箇所）
から局所変形を生じた．

36

Stress, σ , GPa

30

(e)
(d)

1

zz

20

(c)

(b)
0.5

10

(a)
0

(f)

0.02

0.04

0.06

0
0.08

Strain, εFig. 6.1



0

-0.2

-0.4

-0.6

-0.8

Minimum eigenvalue, η (1)min, GPa

1.5

Maximum deviation of cross-sectional centerfrom initial position, r, ÅStress, zz
Maximum deviation of cross-sectional centerfrom initial position, r
Minimum eigenvalue, η (1)min
40

Stress-strain curve，deviation from axis center and change in the minimum of eigenvalue under [1̄010] compression.

Stress, σ , GPa

(e)
30
1

(d)

z
z

20

(c)
0.5

10

(b)
(a)
0

0.02

(f)
0.04

0.06

Strain, εFig. 6.2

0
0.08



1500

1000

500

Number of unstable atoms, Ninst

1.5

Maximum deviation of cross-sectional centerfrom initial position, r, ÅStress, zz
Maximum deviation of cross-sectional centerfrom initial position, r
Number of unstable atoms, Ninst
40

0Stress-strain curve，deviation from axis center and the number of η α(1) <0 atoms under [1̄010] compression.

37

z
y

x

(a)

(b)

(c)

(d)

(e)

(f)

(i) Side view

(a)

(b)

(c)

(d)

(e)

(f)

(ii) Top view (Cross section under periodic boundary in z-axis.)
Fig. 6.3

Snapshots of Mg nanowire at Points (a)∼(f) in Figs. 6.1 and 6.2. Atoms
are colored by the sign of η α(1) , green ; η α(1) > 0, red ; η α(1) < 0.

38

6.2

角部原子の前駆挙動

図 6.4 に η α(1) < 0 原子の対応する固有ベクトルの成分の絶対値の平均の変化を示す．図には湾曲度合いを示すグラフ，塑性変形を生じたひずみと，先の図 6.1∼6.3 の(a)∼(f) 点を合わせて示している．[0001] 圧縮，[1̄21̄0] 圧縮と異なり，湾曲が増加す
るよりずっと前のひずみ 0.01 で顕著な変化を生じている．特に |∆γyz | の変化が大き
く，ひずみ 0.01 を境に急減から急増に転じる．|∆γzx | も 0.01 までは著しく増加し，固有ベクトル成分の中で最も高い値を示している．垂直成分も |∆εyy | の減少，|∆εzz |の増加を生じている．ずれはじめる (c) 点付近では角部の η α(1) < 0 原子のモードに顕著な変化は見られない．|∆γxy | はひずみ 0.04 まで一定で変化が少ない．

Eigenvector,

(b)

0.8

40

0.6

30

0.4

20

0.2

10

0

0.02

0.04

Strain, εzz
Fig. 6.4

0.06

0
0.08

Maximum deviation of cross-sectional centerfrom initial position, r, Å(a)| γyz|| γzx|| γxy|| εxx|| εyy|| εzz|(c)(d) (e)(f)Change in the components of the 1st eigenvector of η α(1) < 0 atoms anddeviation from the axis center under [1̄010] compression. Points (a)∼(f)correspond to those in Figs. 6.1∼6.3.

39

6.3

角部原子の不安定モード

図 6.5 に図 6.3(a) のひずみにおける yz 表面と，図 6.3(a)∼(c) のひずみにおける xz
表面および xy 断面のスナップショットを拡大して示し，角部原子の不安定モードを重ねて示した．本圧縮でも図 4.8 で示したような η α(1) < 0 の原子列のずれは生じていない．(i)xz 表面では裏側のモードが見えないが，(ii)xy 断面を見ると対角頂点のモード（1 と 3，2 と 4）が同じであるとわかる．圧縮前の不安定モードは xy 断面でみると放射状であるが，xz 表面，yz 表面を見ると 1 の原子列のモードはほとんど斜め上を向いている．他の 2 列（2 と 4）においてもモードの方向に定まりはないが斜めを向いているものが多い．これは hcp の結晶方位における隣接面の空きサイトへの方向を示している．(a) → (b) で xz 表面では 1 の原子列は ±x 方向モード，4 の原
子列はある程度方向の定まった斜めのモードになっている．xy 断面の図 (ii)(b) で 1と 3 原子列はモードが 3 方向に限定され，その中でも横方向のモードが多く発生している．2 と 4 の原子列は横方向が増えるものの依然として放射状に発生しており，
そのモードは小さい．湾曲し始める (c) では xy 断面の図 (ii)(c) より x 方向の矢印が大きくなっている．
図 6.6 に (b) 点における不安定モードと単位格子の関係を示した．左の六角柱の中の 1∼4 で示した原子を xy 面に投影したのが左右の模式図で，角部 1∼4 の原子列の単位格子での位置と不安定モードの方向を示した．モードの大きさが同じなのは[1̄21̄0] 圧縮と同様に 1∼4 の原子が等価だからである．

40

2

1

4

1

(a)

4

1

(a)

4

1

(b)

(c)

(i) Side view
2

3 2

3 2

3

1

4 1

4 1

4

(a)

(b)

(c)

(ii) Top view (Cross section under periodic boundary in z-axis.)
Fig. 6.5

Distribution of η α(1) <0 atoms and unstable mode in Fig. 6.3.

2

2

z
1

Fig. 6.6

y

3

3

1
4

x

1

Position in hexagonal prism of the four edge lines in [1̄010] Mg wire.

41領域毎の AES の固有値6.4

図 6.7 は図 4.4 と同様に，図中右の断面で色分けした場所ごとに η α(1) を平均して圧縮時の変化を示したものである．η α(1) はどの場所でも局所変形を生じるひずみ 0.05までは減少しているが，(1̄21̄0) 表面のみ他の場所より減少が小さく湾曲が増加する
少し前のひずみ 0.025 で (0001) 表面と入れ替わる．

Eigenvalue, η (1)ave, GPa

Inn

er
30

10



20
5

Corne
r
10

0

0.02

0.04

0.06

0
0.08

Maximum deviation of cross-sectional centerfrom initial position, r, Å

40

15Conner

Inner

Strain, εzz(a) Average of η α(1) for each zone in (b) (b) Zone in cross section for η α(1) evaluaand deviation from axis center.
tion.

Fig. 6.7

Change in η α(1) averaged by conner, surface and inner under [1̄010] compression.

42

6.5

領域毎の固有ベクトル成分の変化[0001] 圧縮，[1̄21̄0] 圧縮では角部の η α(1) < 0 原子に着目して議論を行っていたが，本圧縮においては座屈挙動との対応が見られなかった．そこで角部の η α(1) < 0 原子を含めたすべての原子に対して η α(1) の固有ベクトルの成分を評価した．図 6.8 に図

6.7(b) の領域によって分けて評価した固有値 η α(1) の，対応する固有ベクトルの成分の絶対値の平均の変化を示す．η α(1) < 0 原子が前駆挙動を示したひずみ 0.01 では
角部の 4 原子列で平均した黒線が顕著な変化を示すが，緑線の (1̄21̄0) 表面の原子も|∆εzz |，|∆γyz | が似た挙動を示している．湾曲が増加するひずみ 0.03 近傍で前駆的な変化を示しているのはやはり |∆εzz | と |∆γyz | で，特に青色の (0001) 表面が |∆γyz |が著しく低下して |∆εzz | が急増している．この (0001) 表面は |∆γzx | でも特徴的な変化を示しており，圧縮により著しく増加している．内側の原子で最も高いのがこの |∆γzx | のせん断モードであるが，これは底面すべりに対応すると考えられる．表面の原子面はバルクよりも収縮しているので圧縮前は低いが，圧縮による横ひずみ
で原子間隔が広がって |∆γzx | 方向の変形モードが大きくなったものと考える．

43

Corner

Inner
Maximum deviation of cross-sectional centerfrom initial position, r

10

0.2

0.02

0.04

0.06

0
0.08

Strain, εzz|

10Eigenvector,

0.2

0.04

0.06

0
0.08

Strain, εzz|

Eigenvector,

20 

0.4

10

0.2

0

0.02

10

0.2

0.06

Strain, εzz

0.04

0
0.08

0.06

0
0.08

Strain, εzz

|
40

0.8

30

0.6

20 

0.4

|

30

0.6

20

0.4

10

0.2

0

0.02

0.04

0.06

Strain, εzz

|

|

0
0.08

Maximum deviation of cross-sectional centerfrom initial position, r, Å

30

0.6

Fig. 6.8

40

Maximum deviation of cross-sectional centerfrom initial position, r, Å

40

0
0.08

|

|

0.8

0.06

0.8

|

0.04

0.04

Maximum deviation of cross-sectional centerfrom initial position, r, Å

20

0.4

0.02

0.02

Strain, εzz

Maximum deviation of cross-sectional centerfrom initial position, r, Å

30

0.6

0

0

|
40

0.02

10

|

0.8

0

20

0.4

0.2

Eigenvector,

0

30

0.6

Eigenvector,

20

0.4

40
0.8

Eigenvector,

Eigenvector,

30

0.6

Maximum deviation of cross-sectional centerfrom initial position, r, ÅMaximum deviation of cross-sectional centerfrom initial position, r, Å

40
0.8



|

Change in the components of the eigenvector for each zone in Fig. 6.7(b)and deviation from axis center under [1̄010] compression.
44

7

結言
分子動力学法により正方形断面の単結晶 Mg ナノワイヤの座屈シミュレーション

を 3 つの結晶方位について行い，ワイヤが湾曲し始める点を原子弾性剛性係数マト
リクス Bijα =∆σiα /∆εj の第一固有値 η α(1) とその固有ベクトルの観点から議論した．
得られた結果を以下に示す．

4 章では [0001] 方向の圧縮について解析を行った．その結果，
1. 湾曲の開始を定量的に評価するため，軸方向の各原子面（角柱断面）の重心位
置 (xi , yi ) の移動量を測定し，その最大値をモニターした．その結果 [0001] 圧
縮ではひずみ 0.02 付近で湾曲し始め，ひずみ 0.03 近傍から急激に増加していることを明らかにした．

2. 圧縮前から正方形断面の四隅の 1 原子列が η α(1) < 0 であった．湾曲し始める
ときに新たに η α(1) < 0 の原子が発生することはなかった．

3. 断面の角部（corner），x 方向表面（(1̄21̄0) 面），y 方向表面（(1̄010) 面），内部（inner）の 4 つの領域に分け，場所ごとに η α(1) を平均して圧縮時の変化を調べたところ，角部以外では固有値が増加（安定化）していた．湾曲開始点近傍で固有値の増加は停止した．

4. 2 の角部の η α(1) < 0 原子の固有ベクトルの成分の絶対値 |∆εxx |, |∆εyy |, · · · , |∆γxy |の平均変化を調べたところ，湾曲が増えはじめる直前で成分の大小関係の急激な入れ替わり＝変形モード変化が湾曲の前駆挙動として確認された．

5. η α(1) < 0 原子の固有ベクトル {∆ε1 , ∆ε2 , · · · , ∆ε6 }={∆εxx , ∆εyy , · · · , ∆γxy } から最大せん断方向をワイヤ側面や断面上の矢印として示し，湾曲時のモードを可視化した．

45

5 章では [1̄21̄0] 方向の圧縮について解析を行った．その結果，
6. この方向の圧縮では，塑性変形（原子構造変化）を生じないままひずみ 0.05をピークとする滑らかな正弦波状の応力ーひずみ曲線を示した．

7. 四隅ではなく突き出た対角の 1 原子列が圧縮前から η α(1) < 0 であり，6 のピー
ク近傍で新たに η α(1) < 0 となる原子が発生した．

8. 湾曲はひずみ 0.03 近傍から始まっていた．このとき η α(1) < 0 の原子数に変化はない．

9. 断面の領域毎に平均した η α(1) はいずれも圧縮によって減少し，湾曲が急増しはじめると極小値を示して増加に転じた．

10. 角部の η α(1) < 0 原子の固有ベクトル成分（絶対値）の変化を調べたところ，[0001] 圧縮と同じようにモード変化が湾曲し始める前駆挙動として生じていた．

6 章では [1̄010] 方向の圧縮について解析を行った．その結果，11. [0001] 圧縮と同じように応力ピーク後，塑性変形発生による応力急減を生じた．
湾曲はひずみ 0.03 近傍から開始していた．

12. 圧縮前から角部の四隅の 1 原子列が η α(1) < 0 であり，湾曲開始時は η α(1) < 0原子数に変化はない．

13. 4 章，5 章と異なり，角部の η α(1) < 0 原子のモード変化は湾曲が増加するずっと前のひずみ 0.01 で生じており，湾曲し始める前駆挙動とは考えられない．

14. 断面領域毎に η α(1) > 0 の原子の固有ベクトルも含めて成分変化を調べたとこ
ろ，湾曲が増加するひずみ 0.03 の直前に (0001) 表面，ついで内側のバルク部の |∆εzz | と |∆γyz | が急変する前駆的な挙動を示していた．

46

参考文献
(1) 屋代, 西村, 冨田, 材料, 54(10), pp.1053-1059, (2005).
(2) 屋代, 西村, 冨田, 材料, 57(2), pp.112-118, (2008).
(3) 屋代，機械の研究, 養賢堂, 60(7), pp.743-753, (2008).
(4) K. Yashiro, Computational Materials Science, 131, pp.220-229, (2017).
(5) K. Yashiro, et al., Computational Materials Science, 183, 109804, (2020).
(6) J. Wang, et al., Physical Review Letters, 71(25), pp.4182-4185, (1993).
(7) J. Wang, et al., Physical Review B, 52(17), pp.12627-12635, (1995).
(8) 太田，屋代，内藤， Mg ナノワイヤの座屈シミュレーションと原子弾性剛性係数による変形モード解析 ，卒業論文，(2020).

(9) C. Y. Wang, et al., International Journal of Solids Structures, 40(15), pp.38933911, (2003).
(10) X. Han, et al., Composites Science and Technology, 68(7-8), pp.1809-1814,(2008).
(11) S. Ajori, and R. Ansari, Current Applied Physics, 14(8), pp.1072-1077, (2014).
(12) M. Goel, et al., Technical article-peer-reviewed, 20, pp.868-881, (2020).
(13) Z. Wang, et al., Physical E: Low-dimensional Systems and Nanostructures,
40(3), pp.561-566, (2008).
(14) Y. Jing, et al., Computational Materials Science, 45(2), pp.321-326, (2009).

47

(15) M. Nishimura, et al., The Japan Society of Mechanical Engineers, 26, 1509,(2013).
(16) M. Nishimura, et al., Key Engineering Materials, 626, pp.234-239, (2015).
(17) M. Nishimura, and N. Takahashi, The Japan Society of Mechanical Engineers,
53, 1113, (2016).
(18) 上田顯, コンピュータシミュレーション, 朝倉書店,（1990）．
(19) 洲之内治男, サイエンスライブラリ–理工系の数学=15, 数値計算, サイエンス社,（1978）．

(20) 香山正憲, 固体材料の電子状態の基礎, 材料学会勉強会資料, (1993)．
(21) X. W. Zhou, et al., Physical Review B, 69, 144113, (2004).
(22) X. W. Zhou, et al., Physical Review B, 69, 035402, (2004).
(23) X. W. Zhou, et al., Acta materialia, 49(19), pp.4005-4015, (2001).
(24) D. C. Wallace, Thermodynamics of Crystals, Wiley, Newyork, (1972)．
(25) K. Yashiro, AIP Advances, 10, 035301, (2020).

48

謝辞
本研究を遂行するにあたり，無知無学な著者に対し多くの丁寧なご指導を賜りました屋代如月教授に心より御礼申し上げます．研究に関して様々な視点から御助言を頂いた内藤圭史助教にも厚く感謝を申し上げます．また，共に研究を進め切磋琢磨し合った同輩の研究室メンバー，これまで数多くの御助言を頂きました先輩方，後輩方に心より感謝いたします．

49