2017.2.13
修士論文審査会	

原子スケールにおける	  異種材界面の応力場・ひずみ場評価	
鹿児島大学大学院 理工学研究科
機械工学専攻 固体力学研究室
修士課程2年 芝健太	

発表の流れ	• 研究背景・目的• 解析理論-異方性弾性論-原子シミュレーション理論• 解析条件・結果-分子静力学解析による応力場・ひずみ場-MS解と弾性解の比較• 結言	

発表の流れ	• 研究背景・目的• 解析理論-異方性弾性論-原子シミュレーション理論• 解析条件・結果-分子静力学解析による応力場・ひずみ場-MS解と弾性解の比較	• 結言	

研究背景・目的	
電子デバイスの微細化や積層化の著しい発展	
Material B	
材料が積層され，多くの異種材界面が存在	

格子定数の違いにより接合界面で応力集中が発生	

Material A	

破壊や剥離の起点となり，電子機器の破壊や故障の原因	
• デバイスの厚みや薄膜技術がナノスケールで制御	  • 積層されている電子デバイスのほとんどは異方性材料	

原子スケールにおける	  異方性異種材界面の応力場・ひずみ場の評価が重要	  研究背景・目的	
材料定数の違いにより応力が集中する原因	  	異種材整合面	
6原子	

6原子	

Material 27格子定数比	 ‥6	大きい弾性ひずみ場（集中荷重が作用）	

6原子	

Material 16原子	研究背景・目的	
材料定数の違いにより応力が集中する原因	  	異種材整合面	
6原子	

6原子	

Material 2

圧縮	  圧縮	  7格子定数比	 ‥6	引張	  引張	  7原子	Material 17原子	研究背景・目的	
材料定数の違いにより応力が集中する原因	  	異種材半整合面	
6原子	

6原子	

Material 27格子定数比	 ‥6	弾性ひずみを緩和	  7原子	Material 1

7原子	

ミスフィット転位（弾性ひずみを緩和するための転位）	

研究背景・目的	
ミスフィット転位とは・・・	
異種材間の格子定数が異なることが原因で生じる	  弾性ひずみを緩和させるため，界面に周期的に発生する転位	

Material 2Material 1

透過型電子顕微鏡（TEM）画像[1]格子定数が異なる異種材界面（半整合面）	
[1] A. Jallipalli, et al., Nanoscale research letters (2009).	

研究背景・目的	
先行研究	
森山らが仮想モデルにおける異方性異種材界面・角部で，• 分子静力学解析より，ミスフィット転位が破壊靭性値に	  及ぼす影響を評価	  • 分子静力学法(MS法）と異方性弾性解の結果を比較し，格子定数差が小さいときは両者が比較的一致していたと，ミスフィット転位の影響を評価 [2]．

目的	
実材料に生じるミスフィット転位の応力場・ひずみ場を弾性論と解析から定量的に評価→積層材料の信頼性向上を図る[2] 森山真樹，ミスフィット転位が存在する異方性異種材界面および界面角部の応力・ひずみ場，京都大学修士論文，(2014)．	

研究背景・目的	
評価方法	
比較	
異方性弾性論[3] による	  応力場・ひずみ場の評価	  原子シミュレーションを用いた解析	  • 先行研究では仮想モデルのみの評価であり，実材料を仮定していない	• 弾性論の見地からの評価は積層材料の信頼性向上に必要• Ge/Siブロックモデルを作成し，原子シミュレーション（分子静力学）から応力・ひずみ場を算出．

[3] C. Hwu, Anisotropic Elastic Plates, Springer Science+Business Media, (2010), pp. 53-113.	

発表の流れ	• 研究背景・目的• 解析理論-異方性弾性論-原子シミュレーション理論• 解析条件・結果-分子静力学解析による応力場・ひずみ場-MS解と弾性解の比較	
• 結言	

解析理論（異方性弾性論）	
ミスフィット転位を含む異種材界面（半整合面）の応力	
σ semicoherent = σ coherent + ∑σ dislocation
半整合面における応力	 整合面における応力	 転位周りの応力の重ね合わせ	整合面(coherent)	

Material 2

整合面における応力	

σ incoherent = σ coherent + ∑σ dislocation
半整合面(semicoherent)	

σ incoherent = σ coherent + ∑σ disloca
Material 1
Material 2

ミスフィット転位周りの応力σ incoherent = σ coherent + ∑σ dislocation σ incoherent = σ coherent + ∑σ dislocation
Material 1

整合面における応力（異方性弾性論）	
整合面における応力σcoherent	 条件	・側面の応力は接合面から生じている 	σ A +σ B = 0

・上下面は自由	σ yA = σ yB = 0

・材料A, Bの一辺の長さは等しくなる 	a(1+ ε A ) = b(1+ ε B )・フックの法則が成立する	
εA =1(σ A − ν Aσ A )
EA

εB =1
(σ B − ν Bσ B )EBE:	  ヤング率	  ν:	  ポアソン比	  S:弾性コンプライアンス	

単位格子で考える	(ax − bx )(a z S zzA + bz S zzB ) − (a z − bz )(ax S xzA + bx S xzB )σ = −σ =(ax S xzA + bx S xzB )(a z S xzA + bz S xzB ) − (ax S xxA + bx S xxB )(a z S zzA + bz S zzB )Ax

Bx

等方性材料と考えると	

σ incoherent = σ coherent + ∑ dislocation

σσ A = −σ B =EA EB (b − a)
(EA + EB )(1− ν )転位周りの応力（異方性弾性論）	
ミスフィット転位における応力σdislocation		境界条件(Boundary Condition)（転位を原点とする極座標系）	

u1 (r,α 0 ) = u 2 (r,α 0 ),φ1 (r,α 0 ) = φ2 (r,α 0 ) u:変位

Φ:応力関数p：集中荷重b:バーガースベクトルh, g:実数ベクトル
Ñ(k) : 基本弾性行列	

φ1 (r,α 0 + π ) − φ2 (r,α 0 − π ) = p̂u1 (r,α 0 + π ) − u 2 (r,α 0 − π ) = b̂,∀ r

σ ij → 0,at infinity	 (k ) (θ ) − N
 (k ) (α )]h +[ N
 (k ) (θ ) − N
 (k ) (α )]g}
u k = ln rh + π {[ N110
220(k)(k)(k)T(k)T (θ ) − N
 (α )]h +[ N (θ ) − N (α )]g}
φ = ln rg + π {[ Nk330110

極座標(r,θ)における応力と応力関数の関係式より	
tθ = φ,rt r = φ,θ / rnT g
π sT (k )σ θθ =,σ rr = −[N 3 (θ )h + N1(k )T (θ )g]
rr
TTs gπnσ rθ ==−[N 3(k ) (θ )h + N1(k )T (θ )g]rr"h%πσ θ 3 = $ ' ,σ r3 = − [N 3(k ) (θ )h + N1(k )T (θ )g]3r# r &3σ rr = −sT φ,θ / r σ rθ = sT φ,r = −nT φ,θ / rincoherent
σ r3 = −(φ,θ )3 / r σ θ 3 = (φ,r )3σ θθ = nT φ,rσ= σ coherent + ∑σ dislocation

転位周りの応力（異方性弾性論）	
ミスフィット転位における応力σdislocation	
−1 T−1T−1 TT	
N
=
TR,N
=
T
=
N,
N
=
RT
R−
Q=N1 Condition) 22
33,境界条件(Boundary"C11 C16 C15 %"C16 C12 C14 %"C66（転位を原点とする極座標系）	
$'$'$Q=CCC,R=C
CC,T=u1 (r,α 0 ) = u 2 (r,α 0 ),φ1$(r,16α 0 ) =66 φ2 (r,56α' 0 ) u:変位
2646 '$ 66$C26
Φ:応力関数$'$#C56 C25 C45 '&$#C46CCC#&155655
φ (r,α + π ) − φ (r,α − π ) = p̂
102p：集中荷重
b:バーガースベクトルh, g:実数ベクトル
Ñ(k) : 基本弾性行列	

0u1 (r,α 0 + π ) − u 2 (r,α 0 − π ) = b̂,∀ r

σ ij → 0,€

at infinity	 (k ) (θ ) − N
 (k ) (α )]h +[ N
 (k ) (θ ) − N
 (k ) (α )]g}
u k = ln rh + π {[ N110220(k)(k)(k)T(k)T (θ ) − N
 (α )]h +[ N (θ ) − N
 (α )]g}
φ = ln rg + π {[ Nk330110

極座標(r,θ)における応力と応力関数の関係式より	
tθ = φ,rt r = φ,θ / r

σC24C46 %'C24 'C44 '&nT g
π sT (k )σ θθ =,σ rr = −[N 3 (θ )h + N1(k )T (θ )g]
rr
TTs gπnσ rθ ==−[N 3(k ) (θ )h + N1(k )T (θ )g]
rr"h%πσ θ 3 = $ ' ,σ r3 = − [N 3(k ) (θ )h + N1(k )T (θ )g]3r# r &3σ rr = −sT φ,θ / r σ rθ = sT φ,r = −nT φ,θ / r
incoherent
σ r3 = −(φ,θ )3 / r σ θ 3 = (φ,r )3σ θθ = nT φ,rC26
C22= σ coherent + ∑σ dislocation

発表の流れ	• 研究背景・目的• 解析理論-異方性弾性論-原子シミュレーション理論• 解析条件・結果-分子静力学解析による応力場・ひずみ場-MS解と弾性解の比較	
• 結言	

分子静力学法（原子シミュレーション理論）	
• 原子間のポテンシャルによって相互作用する質点で近似し，構造緩和を行うことで系全体のポテンシャルエネルギーを最小化することで系の特性や物性を求めることができる．
• 0Kでの分子動力学法に対応し，静的に原子構造を取り扱う．
→理想的な系での計算を行うことができる．	

ポテンシャル	

ポテンシャル	fiα = −原子間距離	∂Φ∂riαfiα：原子αのi方向成分の力
Φ：全ポテンシャルエネルギー
riα：原子間距離	

原子間距離	

原子間ポテンシャル（原子シミュレーション理論）	
原子間ポテンシャル → Stillinger-Weberポテンシャル[4]を使用	
二体間ポテンシャルに三体間ポテンシャルを付加した形で表現されるポテン
シャルエネルギーで，ダイヤモンド構造のパラメータにフィッテングしている．	

全ポテンシャルエネルギーΦ	Φ = ∑∑ v2 (i, j) + ∑∑∑ v3 (i, j, k)
ij>i

ij>i k> jv2：二体間のポテンシャルエネルギーv3：三体間のポテンシャルエネルギー
rij,rik：原子iと原子j, kとの距離
θijk：原子i-jと原子i-kの結合角度k
rik

θijk

irij

[4] Frank H. Stillinger and Thomas A. Weber, Computer simulation of local order in condensed
phases of silicon, Physical Review B, Vol. 31, Num. 8, (1985), 5262.	

j原子間ポテンシャル（原子シミュレーション理論）	
l 二体間ポテンシャル	
v2 (rij ) = ε f2 (rij / σ )! A(Br − p − r −q )exp[(r − a)−1 ] , r < af 2 (r) = "#0 , r ≥ al 三体間ポテンシャル	

v3 (ri , rj , rk ) = ε f3 (ri / σ , rj / σ , rk / σ )f 3 (ri ,r j ,rk ) = h(rij ,rik ,θ jik ) + h(rji ,rjk ,θ ijk ) + h(rki ,rkj ,θ ikj )1 2h(rij ,rik ,θ jik ) = λ exp[γ (rij − a) + γ (rik − a) ]× (cosθ jik + )3−1−1ε（eV）：結合エネルギーの決定するパラメータσ（Å）：原子間距離を決定するパラメータ
完全結晶のみ適用可能なため，
f2, f3：ポテンシャル関数
cosθ=-1/3を使用（θ=109.5°）	
A, B, p, q, a, λ, γ：ポテンシャルパラメータ	

力の応力変換式	
fiα = −∂Φ∂riα

Φ = ∑∑ v2 (i, j) + ∑∑∑ v3 (i, j, k)
ij>i

ij>i k> j

力fを以下の式で応力に変換	

1 ∂(K + φ )1∂φ ∂r αβα α ασ ij = − () = − (∑ m vi v j + ∑ αβ)V∂εij
V α∂εij
α <β ∂r=−1(∑ mα viα vαj + ∑ (riα − ri β ) f jαβ )V αα <βN αpσ ijα = −11α α αα αββ βα(mvv+(rf+r∑i jiji f j ))
αV
2 n=1

ひずみ変換	

K：運動エネルギー m：質量 v：速度 Vα：原子αの体積
Npα：原子αの配位数 σij：系全体の応力 σijα：原子αの応力 	

C11(GPa)	 C12(GPa)	 C44(GPa)	

応力のひずみ変換式	

Si [5]	

166	

63.9	79.6	

Ge [5]	

129	48.3	66.8	

分子静力学解析から得られた応力場をひずみ場に変換する方法としてHooke’s法則を使用した．	
σ：応力σ ij = Cijklε kl

C：弾性スティフネス
ε：ひずみ	εij = Sijklσ kl

S：弾性コンプライアンス

これをεについて解くとSiやGeのような立方晶を持つ材料は，S11，S12，S44の弾性スティフネスから表される	!ε xx $ ! S11# & ##ε yy & ##ε & ## zz & = ##ε xy & ## & ##ε yz & ##"ε zx &% #"S12S1200

S11S1200

S1100

S440Sym.

S440 $!σ xx $&# &0 &#σ yy &&0 &#σ zz &# &0 &#σ xy &&# &0 &#σ yz &S44 &%#"σ zx &%[5] R. Srinivasan, Lattice theory of third order elastic constants of Germanium and Silicon, J. Phys. Chem.
Solids Pergamon Press, Vol. 28, (1967), pp. 2385-2369.		

発表の流れ	• 研究背景・目的• 解析理論-異方性弾性論-原子シミュレーション理論• 解析条件・結果-分子静力学解析による応力場・ひずみ場-MS解と弾性解の比較	
• 結言	

解析条件・モデル	
y [010]aGe×48aGe×30aSi×30z [001]

Lattice size and lattice size gap.	

aGe×48Ge

lattice size
Si (Å)	

lattice size
Ge (Å)	

lattice sizegap (%)	

5.43	

5.65	

4.1	

Model size.	

Si

aSi,Ge (nm)	 aSi+aGe (nm)	 aSi,Ge (nm)	

x [100]aSi×50ε (eV)	 σ (Å)	27.16	

aSi×5033.27	

27.16	SW Potential parameters for Si and Ge[6].	
λ	p	q	γ	Si	2.17	 2.095	 21.0	 7.050	 0.602	 4.0	

0.0	

1.20	 1.80	 -0.333	

Ge	1.93	 2.181	 31.0	 7.050	 0.602	 4.0	

0.0	

1.20	 1.80	 -0.333	

0	1.2	

Si-Ge	 2.05	 2.138	 26	A	0	B	

0	0	a	cosθ	

1.8	 -0.333	

[6] Zi Jian, Modification of Stillinger-Weber potentials for Si and Ge, Physical Review B, Vol. 41, Num. 18,(1990), 12915.	

解析条件・モデル	
計算手順	
ブロックモデルを作成
x, z方向に周期境界条件を適用	

Boundary Ge	atoms	
Si	
yz

x

y方向の上下2単位格子分を固定	

構造緩和計算→原子が最安定位置になるまで行う
構造緩和計算	
LAMMPSを用いて共役勾配法 (CG法) で行う．	
Diamond structure	

②①

解析結果（MS解）	

Ge③Siyzx10GPa10.01①	②	

①	②	-7.638③	③	σxx	 界面から1原子上下の状態	 σyy	-10GPa-10.29

• ミスフィット転位の発生による応力集中（Ge側で引張力，Si側で圧縮②①力）が発生．
解析結果	
Ge• 転位近傍以外の界面では，ミスフィット転位に作用した応力を緩和するためにGe側で圧縮力，Si側で引張力が生じていた．
③• 界面上ではクロス上に転位線が生じていた．
Siyzx10GPa10.01①	②	

①	②	-7.638③	③	σxx	 界面から1原子上下の状態	 σyy	-10GPa-10.29

解析結果（MS解による応力場）	
y [010]GeSi

σxx	
σyy	σxx	σyy	②：[101]①：[100]x [100]z [001]① [100]の応力分布（1原子上下）② [101]の応力分布（1原子上下）	

σxx	σyy	① [100]	 界面から1原子上（Ge側）	
② [101]	

σxx	
σyy	

σxx	
σyy	

① [100]	 界面から1原子下（Si側）	 ② [101]	
• ミスフィット転位の影響による応力振幅が発生している• ミスフィット転位近傍では±7GPaほど，転位線上には±5GPaほどの応力が作用している• 分子静力学のような理想的な系での計算では，ミスフィット転位が約14nm毎に発生する

解析結果（MS解によるひずみ場）	
y [010]

GeSi

εxx	
εyy	
②：[101]①：[100]x [100]z [001]εxx	
εyy	

① [100]の応力分布（1原子上下）② [101]の応力分布（1原子上下）	

εxx	
εyy	

① [100]	 界面から1原子上（Ge側）	
② [101]	
εxx	
εxx	
εyy	
εyy	

① [100]	 界面から1原子下（Si側）	 ② [101]	

• ミスフィット転位近傍では±3~4%ほど，転位線上には±2~2.5%ほどのひずみが作用している• SiよりもGeのほうがひずみ量が大きい（ヤング率の影響によるもとだと考えられる）．

解析結果（MS解によるひずみ場）	
y [010]

GeSi

εxx	
εyy	
②：[101]①：[100]x [100]z [001]

GeSi
εxx：3.3% -2.7%
εyy：3.7% -2.5%
εxx	
εyy	

① [100]の応力分布（1原子上下）② [101]の応力分布（1原子上下）	

εxx	
εyy	

① [100]	 界面から1原子上（Ge側）	
② [101]	
εxx	
εyy	

εxx	
εyy	

① [100]	 界面から1原子下（Si側）	 ② [101]	

• ミスフィット転位近傍では±3~4%ほど，転位線上には±2~2.5%ほどのひずみが作用している• SiよりもGeのほうがひずみ量が大きい（ヤング率の影響によるもとだと考えられる）．

発表の流れ	• 研究背景・目的• 解析理論-異方性弾性論-原子シミュレーション理論• 解析条件・結果-分子静力学解析による応力場・ひずみ場-MS解と弾性解の比較	
• 結言	

弾性解算出方法	
ミスフィット転位を含む異種材界面（半整合面）の応力	
σ semicoherent = σ coherent + ∑σ dislocation
半整合面における応力	 整合面における応力	 転位周りの応力の重ね合わせ	
lattice size（Å）	 C11(GPa)	 C12(GPa)	 C44(GPa)	

z dislocation lines

ll

Si	

5.43	

166	

63.9	

79.6	

Ge	

5.65	

129	

48.3	

66.8	

ax=az=5.65（Å） bx=bz=5.43（Å）	
(ax − bx )(a z S zzA + bz S zzB ) − (a z − bz )(ax S xzA + bx S xzB )σ = −σ =n=3
(ax S xzA + bx S xzB )(a z S xzA + bz S xzB ) − (ax S xxA + bx S xxB )(a z S zzA + bz S zzB )n=3n=2n=2n=1σcoherentxn=1n=1n=1= 3.14 GPa	n=2n=2n=3 n=3AxBx

格子定数が小さい材料の最小原子間距離をバーガースベク
トルとし，1000個の転位を重ね合わせることで解析を実施．	

Analysis model(bi-material interface)解析結果（MS解と弾性解の比較）	

①	15GPa

15GPa

10.01

10.01

②	③	-7.638③	-15GPa-10.29

分子静力学	

分子静力学	
弾性解	σxx	 弾性解	
σ	yy
界面から1原子上の状態（Ge側）	

-7.638-15GPa-10.29

解析結果（MS解と弾性解の比較）	
• 分子静力学解析と弾性解の結果が定性的に一致している．

①	15GPa

15GPa

10.01

10.01

②	③	-7.638③	-15GPa-10.29

分子静力学	

分子静力学	
弾性解	σxx	 弾性解	
σ	yy
界面から1原子上の状態（Ge側）	

-7.638-15GPa-10.29

解析結果（MS解と弾性解の比較）	

y [010]Ge

plot line

Six [100]σxx	 ①界面から1原子上（Ge側）	 σyy	z [001]① 界面上部の応力分布（1原子上）② 界面下部の応力分布（1原子下）	
分子静力学解析	

σxx	 ②界面から1原子下（Si側）	 σyy	

弾性解	  ・σxx，σyyともに転位芯近傍以外で応力分布がほぼ一致している→誤差が1GPa程度と比較的精度よく応力場を表現できている
分子静力学では転位芯が大きさを持つが，弾性解では転位芯で応力が発散してしまうことにより，転位芯の大きさを考慮できていないことが原因である	

発表の流れ	• 研究背景・目的• 解析理論-異方性弾性論-原子シミュレーション理論• 解析条件・結果-分子静力学解析による応力場・ひずみ場-MS解と弾性解の比較	
• 結言	

結言	
分子静力学解析を用いて，格子定数差4.1%を持つGe/Si界面に発生するミスフィット転位近傍の応力場・ひずみ場を評価した．
以下に得た知見を示す．
• 分子静力学のような理想的な系での計算では，ミスフィット転位による応力集中が発生し，Ge側で引張力，Si側で圧縮力が作用していた．
• ミスフィット転位近傍には，応力が約±7GPa，ひずみに換算すると約±3~4%作用していた．
• MS解と弾性解を比較し，転位芯近傍以外で誤差が1GPa程度であり，応力場が比較的一致していた．