き裂の進展・発生に関する局所格子不安定性解析

修 士 論 文

(

題

目

)

き裂の進展・発生に関する

局所格子不安定性解析

平成

29

年度

岐阜大学大学院

工学研究科 博士前期課程

機械システム工学専攻

目 次

第 1 章 緒言 . . . 1 第 2 章 分子動力学法 . . . 4 2.1 分子動力学法の概要 . . . . 4 2.1.1 EAM ポテンシャル . . . . 5 2.1.2 速度スケーリング法 . . . . 6 2.1.3 弾性剛性係数と格子不安定性 . . . . 8 2.1.4 原子応力 . . . . 8 2.1.5 原子弾性係数 . . . 10 2.1.6 原子弾性剛性係数 . . . 10 2.1.7 局所格子不安定性における原子の変形モード . . . 11 第 3 章 相互作用遮蔽シートを用いたき裂進展に関する検討 . . . 13 3.1 シミュレーション条件 . . . 13 3.2 シミュレーション結果および考察 . . . 15 3.2.1 応力ひずみ曲線 . . . 15 3.2.2 原子応力の分布 . . . 17 3.2.3 不安定解析と AES の固有値に関する検討 . . . 21 第 4 章 ねじり粒界モデルに関する検討 . . . 30 4.1 シミュレーション条件 . . . 30 4.2 シミュレーション結果および考察 . . . 32

4.2.1 Σ5 ねじり粒界の結果 . . . 32 4.2.2 Σ13 ねじり粒界の結果 . . . 38 4.2.3 Σ17 ねじり粒界の結果 . . . 44 第 5 章 円孔の応力集中部からのき裂発生に関する検討 . . . 49 5.1 シミュレーション条件 . . . 49 5.2 シミュレーション結果および考察 . . . 51 5.2.1 応力ひずみ線図および最小固有値の変化 . . . 51 5.2.2 固有ベクトルによる検討 . . . 55 第 6 章 結言 . . . 58 参考文献 . . . 60  

第

1

章 緒言

車体の振動の解析に用いられる有限要素法や車体周りの気流の流れを解析する CFD，新材料の開発には第一原理計算など，工業に使われる計算科学技術の発展は 著しい．この科学技術計算を支えるのが CPU や計算高速化のための GPGPU に用 いられる GPU などのハードウェアであるがその高速化・高機能化はナノレベルの 微細加工技術によって実現されている．最小配線の寸法は数ナノレベルで微細化さ れており，有限要素法等に代表される連続体近似に基づくシミュレーションモデル の適用ができない領域に達している．  ナノレベルの材料挙動をシミュレーションする手法には粒子法 (MPS) や分子動力 学法 (MD) が挙げられる．粒子法は連続体を有限の粒子として離散化する方法であ るが，物質の挙動を詳しくに理解するには物質の最小構成単位である原子を考える 必要がある．分子動力学法により材料の変形挙動を原子レベルから検討しようとす る試みは古くから行われており，1982 年には弾性体とリンクさせた解析が行われて いる(13)_{．また 2003 年には松本ら}(3)_{が，アモルファス金属のき裂の進展挙動につい} て連続体モデルと MD 法による原子モデルと比較した検討を行っている．こういっ たマルチスケールな解析だけでなく，連続体近似を用いずに MD のみでナノレベル での破壊へ適応するプローチも行われており，須賀ら(1)_{は界面き裂に MD 法を適用} した研究を 1991 年に発表している．破壊力学が適用困難な界面破壊を想定して，負 荷モードの一致などから，MD 法が界面き裂の挙動解析に有効であることを示して いる．1997 年には中谷ら(2)_{が，bcc 結晶のき裂先端からの転位挙動について検討し} ている．Po ら(25)_{はニッケル単結晶の片側き裂進展に対して MD シミュレーション} を行い，Ni における (111) すべり挙動の観察や，き裂進展方向の温度により変化す

ることや臨界応力の結晶方位依存性について報告している．また都留ら(29)_{は Al に}

対して 1 億原子の大規模なせん断 MD シミュレーションを行い，塑性変形の開始が 結晶粒内の転位源の張り出しによるものであり，転位理論と MD で得られた臨界応 力が近い値をとることを報告している．こういった MD は Fe, Si, Al, Ni(30)_などの

汎用的な金属以外にも行われており，Qing ら(26)_{は超硬材料 WC 対してナノインデ} ンテーションを行い，(0001)WC が押し込みにより_{{1¯101} 面に転位が射出され転位} ループすることで押し込み反力が低下することを報告している．久保田ら(30)_は純 鉄の多結晶体に対して MD を適用し，変形速度が遅い場合はき裂の進展に対する抵 抗が粒径が小さくなるほど下がることを報告している． 以上のように，変形・破壊現象を原子レベルから検討する試みは多数なされてきた が，問題となるのは観察した事象に対する定量的な考察である．原子間ポテンシャ ルを用いた動的なシミュレーションは，多かれ少なかれ実際の材料の特性から外れ た「モデル材料」にすぎない．原子間力を電子状態から精密に評価する第一原理計 算分子動力学法では扱える原子数はスーパーコンピュータを用いても 1 万個程度で あり，さらに熱揺動の効果は時間の問題から考慮できない．これまでの研究の多く は，原子シミュレーションの結果を応力拡大係数などの連続体理論と比較した議論 を行っているが，連続体近似の有効性を追認するだけになりがちである．したがっ て連続体近似が適用できない事象を観察した時，定量的な議論をするより基準がな い．  我々のグループでは，原子シミュレーションにおける定量的な評価基準として局所 格子不安定性に着目した検討を行ってきた(7)–(11)_{．非線形弾性における応力ーひずみ} 関係を表す弾性剛性係数 Bij = (∆σi/∆εj) を基に，原子弾性剛性係数 (Atomic Elastic Stiffiness;AES)Bα ij = (∆σiα/∆εj) を定義し，その正値性から局所の変形開始を議論す る．局所不安定については，北村, 梅野ら(4), (5)_{が，原子構造体の全自由度 (N 原子な} ら全体の並進・回転自由度を覗いた M = 3N−6) について系のエネルギーの 2 次導関

接解析することは難しい．一方，Bα

ijは 6×6 マトリックスなので，それを全ての原子

について計算するのは MD における力計算とほぼ変わらず 100 万原子でも実行可能で ある．これまで AES によるき裂への適用は Si, Fe, Mg について行われてきた(7)–(11)_．

固有方程式 Bα ij∆εj = ηα∆εiの解 ηα(k)(k = 1∼ 6) の第一固有値 ηα(1)が負の原子の固 有ベクトル_{{εi} から変形モードをひずみ空間 {∆εxx}, ∆εyy, ∆εzz, ∆εyz, ∆εzx, ∆εxy} で示し，さらにはこのひずみ成分 ∆εijからなるひずみテンソルの主軸から変形モー ドを明らかにすることなどを試みている．本研究では，特にき裂の発生に着目した 検討を bcc-Fe について行う．3 章では原子間力を人為的に遮へいすることでき裂開 口させる理想き裂のシミュレーションを行った．4 章ではき裂発生源として CSL ね じり粒界の引張りシミュレーションを行った．5 章では，応力集中部からのき裂発 生シミュレーションとして円孔の引張りシミュレーションを行った．

第

2

章 分子動力学法

2.1

分子動力学法の概要

分子動力学法 (molecular dynamics method，略して MD 法) は，系を構成する各 粒子についてニュートンの運動方程式 mαd 2_rα dt2 = F α _(2.1.1) をたて，これを数値積分することにより粒子の軌跡を求める方法である(14)_．ここ で，mα_，rα_{はそれぞれ原子 α の質量および位置ベクトルである．原子 α に作用す} る力 Fα_{は，系のポテンシャルエネルギー E} totの各位置における空間勾配として次 式により求められる． Fα =−∂Etot ∂rα (2.1.2) 式 (2.1.1) の数値積分には，Verlet の方法，予測子–修正子法等がよく用いられる(15)_． 本研究では，以下に示す Verlet の方法を用いた．時刻 t + ∆t と t− ∆t での粒子 i の 位置ベクトル rα_{(t± ∆t) を Taylor 展開すると} rα(t + ∆t) = rα(t) + ∆tdr α_(t) dt + (∆t)2 2 d2_rα_(t) dt2 + O ( (∆t)3 ) (2.1.3) rα(t− ∆t) = rα(t)− ∆tdr α_(t) dt + (∆t)2 2 d2_rα_(t) dt2 + O ( (∆t)3) (2.1.4) となる．ここで，vα_{を時刻 t における原子 α の速度とすると，} drα = vα(t) (2.1.5)

であり，式 (2.1.1) と式 (2.1.5) を式 (2.1.3) と式 (2.1.4) に代入すると， rα(t + ∆t) = rα(t) + ∆tvα(t) + (∆t) 2 2 Fα(t) mα + O ( (∆t)3) (2.1.6) rα(t− ∆t) = rα(t)− ∆tvα(t) +(∆t) 2 2 Fα(t) mα + O ( (∆t)3) (2.1.7) となる．両式の和と差をとると， rα(t + ∆t) + rα(t− ∆t) = 2rα(t) + (∆t)2F α (t) mα + O ( (∆t)4) (2.1.8) rα(t + ∆t)− rα(t− ∆t) = 2∆tvα(t) + O((∆t)3) (2.1.9) が得られる．∆t3_{以上の高次項は無視できるとすると，時刻 t + ∆t での位置ベクト} ルと t での速度は rα(t + ∆t) = 2rα(t)− rα(t− ∆t) + (∆t)2 F α (t) mα (2.1.10) vα(t) = 1 2∆t{r α (t + ∆t)− rα(t− ∆t)} (2.1.11) と求められる．t + ∆t での座標を求めるには 2 つの時刻 t と t− ∆t での座標が必要 である．初期の計算 (t = 0) では，t = ∆t での座標 rα_{(∆t) は式 (2.1.6) と初速度か} ら得ることができる．

2.1.1

EAM

ポテンシャル

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

EAM(Embedded atom method；原子埋め込み法) ポテンシャルは金属中の多体効果 を良好に再現することから広く用いられている．密度汎関数理論に基づき，まず金 属材料における系のポテンシャルエネルギー Etotは原子を価電子雲中に埋め込むエ ネルギーと原子間の 2 体間相互作用の和で与えられるとする．さらに，埋め込みエ ネルギーは埋め込む位置の電子密度にのみ依存すると仮定することによって，系全 体のエネルギーは次式のように表わされる． Etot= N ∑ α F ( ¯ρα) + 1 2 N ∑ α N ∑ β(̸=α) ϕ(rαβ) (2.1.12) ここで，¯ραは原子 α の位置における多体効果を考慮する密度を表し，F (¯ρα) は密度 ¯ ρα_{の位置に原子を埋め込むエネルギー，ϕ(r}αβ_{) は距離 r}αβ_{離れた原子 α と β のクー} ロン相互作用である．密度 ¯ρα_{は周囲の原子 β からの寄与 ρ(r}αβ_{) の重ね合わせで与} えられると仮定し ¯ ρα = neighbor_∑ β(̸=α) ρ(rαβ) (2.1.13) で評価する．  本研究では，鉄のシミュレーションで広く用いられている FS(Finnis Sinclar) ポテ ンシャル(17)_{を用いた．FS ポテンシャルでは，} Etot= ∑ α Eα =∑ α { 1 2 ∑ β̸=α ϕ(rαβ)− AF (¯ρα) } (2.1.14) で系のエネルギーを評価する．ここで ¯ ρα = ∑ β̸=α ρ(rαβ) (ρ(rαβ) = (r− d)2+ β(r− d)3/d) (2.1.15) であり，d, β, A はパラメータである．また埋め込みエネルギー F (ρ) は F (ρ) = √ρ と単純な関数形であるのが特徴である．d, β, A の値を表 2.1.1 に示す．

2.1.2

速度スケーリング法

Table 2.1.1 Potential parameters of FS potential. A[eV] d [nm] β 1.828905 0.3569745 1.8 速度スケーリング法は，統計熱力学より導かれる次の原子の運動エネルギーと温度 の関係を用いて，以下のように制御する． 1 2m α_vα iv α i = 3 2kBT (2.1.16) ここで mα_{は粒子 α の質量，v}α i は 温度 T での粒子 α の速度，kBは Boltzmann 定 数で kB=1.38× 10−23[J/K] である．目標の温度 T0 における原子 α の速度を viα0 と おくと vα i0 は式 (2.1.17) のように表される． v_iα 0 = ( 3kBT0 mα )0.5 (2.1.17) 同様に，温度 T の時の原子 α の速度は式 (2.1.18) のように表される． v_iα = ( 3kBT mα )0.5 (2.1.18) 式 (2.1.17) と式 (2.1.18) より以下の式が得られる． vα i0 vα i = ( T0 T )0.5 (2.1.19) つまり，系の温度を T から T0 にするには，式 (2.1.19) の右辺を現在の速度に掛け てやればよい．ただ，これだけでは原子配置に反映されないので，Verlet 法におけ る座標更新の式 (2.1.10) を以下のように置き換える必要がある． rα(t + ∆t) = rα(t) + √ T0 T ( rα(t)− rα(t− ∆t) + (∆t)2 F α (t) mα ) (2.1.20) 平衡状態では，Nose-Hoover 法(20)_{によって得られるカノニカルアンサンブルに一致} することが示されている．

2.1.3

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

応力 σij および弾性係数 Cijkl は，等温過程では σij = 1 V ( ∂F ∂ηij ) , Cijkl = 1 V ( ∂2_F ∂ηij∂ηkl ) (2.1.21) と定義される(21)_{．ここで，F は Helmholtz の自由エネルギー (断熱過程では内部エ} ネルギー U)，V は結晶の体積，ηijは平衡状態 (無負荷とは限らない) からの仮想的 な微小ひずみである．一方，無負荷平衡状態を基準とするひずみ εij と応力 σijの関 係は，2 つの平衡状態間の変形を考えて導出される次の弾性剛性係数によって表さ れる(21)_． Bijkl ≡ (_∆σ ij ∆εkl ) = Cijkl+ (σilδjk + σjlδik+ σikδjl+ σjkδil− 2σijδkl)/2 (2.1.22) ここで，δij はクロネッカーのデルタである．Bijklは非線形弾性における応力-ひず み関係の勾配を表すことから，Wang, Yip らは，ひずみの対称性を考慮したテンソ ル Bsym

ijkl ≡ (Bijkl+ Blkji)/2 の正値性によって結晶の安定性を評価することを提案し

ている(22), (23)_．B ijklの最小固有値が負であれば，対応するひずみに対して応力が負 になる自発的変形モードが存在することになり，不安定と判別される．

2.1.4

原子応力

局所の安定性を評価するための原子弾性剛性係数 Bα ijklの算出に必要な原子応力 σα ijならびに原子弾性係数 Cijklα は，各原子周りの微小ひずみに対するポテンシャル エネルギーの 1 次，2 次変化量として導出される．簡単のため，結晶の内部エネル ギー U が系全体のポテンシャルエネルギー Etot に等しいとする．このとき，応力 は平衡状態からの微小ひずみ η に対するポテンシャルエネルギーの単位体積当たり

ここで，V は平衡状態における系の体積であり，下付添字のローマ文字はテンソル のデカルト座標成分を表す．(2.1.23) 式の微分を求めるため，平衡状態からの仮想 的な均一変形を考える．結晶内の原子 α の位置ベクトルは仮想変形のヤコビ行列 J によって rα = J ¯rα (2.1.24) と変化する．ここで，「¯」は仮想ひずみによる変形前の値を示す．これより，原子 α と 原子 β の間の距離 rαβ _には (rαβ)2 = ¯rαβ_i Gijr¯ αβ j (2.1.25) なる関係が成立する．ただし，Gij = JkiJkj である．仮想変形の Lagrange ひずみテ ンソル ηij は ηij = 1 2 [ Gij − δij ] (2.1.26) であり，その微小量 dηij = 1 2dGij (2.1.27) と式 (2.1.25) の関係から次の関係が得られる． ∂rαβ ∂ηij = r¯ αβ i r¯ αβ j rαβ (2.1.28) 以上の関係を用いると FS ポテンシャルにおける応力は次式で評価される σij = 1 V ∂Etot ∂ηij = 1 V (₁ 2 N ∑ α N ∑ β(̸=α) ∂rαβ ∂ηij ∂Etot ∂rαβ ) = 1 V N ∑ α N ∑ β(̸=α) (₁ 2ϕ ′_(rαβ )− AF′( ¯ρα)ρ′(rαβ) )_rαβ i r αβ j rαβ ここで，Ω = V /N とすれば，各原子位置における原子応力を σα_ij = 1 Ω N ∑ β̸=α (₁ 2ϕ ′_(rαβ_{)− AF}′_{( ¯ρ}α_)ρ′_(rαβ₎ )_rαβ i r αβ j rαβ (2.1.29) と定義すると，系の応力は σij = 1 N N ∑ α σ_ijα (2.1.30)

2.1.5

原子弾性係数

弾性係数も応力と同様に U ≈ Etot の場合には Cijkl = 1 V ∂2Etot ∂ηij∂ηkl (2.1.31) であるので，平衡状態からの仮想均一変形を考えると FS ポテンシャルにおける弾 性係数は以下のようになる． Cijkl = 1 2V N ∑ α N ∑ β(̸=α) ∂rαβ ∂ηkl ∂ ∂rαβ (_∑N α N ∑ β(̸=α) {₁ 2ϕ ′_(rαβ_{)− AF}′_{( ¯ρ}α_)ρ′_(rαβ₎ }_rαβ i r αβ j rαβ ) = 1 V [ 1 2 N ∑ α N ∑ β̸=α { ϕ′′(rαβ)−ϕ ′_(rαβ₎ rαβ }_rαβ i r αβ j r αβ k r αβ l (rαβ₎2 − A∑N α N ∑ β̸=α F′( ¯ρα) { ρ′′(rαβ)− ρ ′_(rαβ₎ rαβ }_rαβ i r αβ j r αβ k r αβ l (rαβ₎2 − A∑N α F′′( ¯ρα){ ∑ β̸=α ρ′(rαβ)r αβ i r αβ j rαβ }{ ∑ β̸=α ρ′(rαβ)r αβ k r αβ l rαβ }] (2.1.32) 応力と同様に，各原子位置における原子弾性係数を以下のように定義する． C_ijklα = 1 Ω [ 1 2 N ∑ β̸=α { ϕ′′(rαβ)−ϕ ′_(rαβ₎ rαβ }_rαβ i r αβ j r αβ k r αβ l (rαβ₎2 − A∑N β̸=α F′( ¯ρα) { ρ′′(rαβ)−ρ ′_(rαβ₎ rαβ }_rαβ i r αβ j r αβ k r αβ l (rαβ₎2 − AF′′_{( ¯ρ}α₎{ ∑ β̸=α ρ′(rαβ)r αβ i r αβ j rαβ }{ ∑ β̸=α ρ′(rαβ)r αβ k r αβ l rαβ }] (2.1.33) これより，系の弾性係数は Cijkl = 1 N N ∑ α C_ijklα (2.1.34) のように原子弾性係数の平均となる．

2.1.6

原子弾性剛性係数

Stiff-Wang らによる提案(22), (23)_{に従い，Voigt 対称性をもたせた B}α sym

ijkl ≡ (Bαijkl+Blkjiα )/2

を用いて安定性評価を行う．なお，ここでの i∼l は直交座標の指標 (1,2,3) もしくは (x，y，z) であり，B_ijklα は 4 階のテンソルであるが，Voigt 対称性を持たせた場合， 弾性係数テンソルと同じく独立な成分は 21 個となり，xx，yy，zz，yz，zx，xy を 1-6 とする Voigt 表記を用いれば 6× 6 のマトリックスとして Bα

ij と表すことができる．

その場合の式 (2.1.35) の対応を表 2.1.2 に示す．

Table 2.1.2 Atomic elastic stiffines matrix in Voigt notation.

Cα 11+ σα1 C12α − σα 1 + σα2 2 C α 13− σα 1 + σ3α 2 C α 14− σα 4 2 C α 15+ σα 5 2 C α 16+ σα 6 2 Cα 22+ σα2 C23α − σα 2 + σ3α 2 C α 24+ σα 4 2 C α 25− σα 5 2 C α 26+ σα 6 2 C33α + σ3α C34α + σα 4 2 C α 35+ σα 5 2 C α 36− σα 6 2 Cα 44+ σ₂α+ σα₃ 2 C α 45+ σα₆ 2 C α 46+ σ₅α 2 Cα 55+ σα 3 + σ1α 2 C α 56+ σα 4 2 sym. Cα 66+ σα 1 + σ2α 2

2.1.7

局所格子不安定性における原子の変形モード

AES の各固有値に対する固有ベクトル{∆εi} は ∆σ_iα = B_ijα∆εj == ηα(n)∆εi = ηα(n)             ∆ε1 ∆ε2 ∆ε3 ∆ε4 ∆ε5 ∆ε6             = ηα(n)             ∆εxx ∆εyy ∆εzz ∆εyz ∆εxz ∆εxy             と書ける．上記のように固有ベクトル_{{∆εi} は 6 次元のひずみ空間におけるベクトル} であり，変形モードを想定するのは困難である．そのため，_{{∆εi} の成分 ∆εxx}, ∆εyy, ∆εzz, ∆εyz, ∆εzx, ∆εxyを成分とする以下のひずみテンソル [∆εij] の主軸を議論する．

[∆εij] =     ∆εxx ∆εxy ∆εxz ∆εyy ∆εyz sym. ∆εzz     ηα(1) _{< 0 の原子の固有ベクトルの x} 1, x2, x3について，[∆εij] の最も大きい固有値に 対応した x3を原子の主変形モードとして議論した．Bijα, [∆εij] の固有値・固有ベク トル解析には数値計算ライブラリ LAPACK を使用した．

第

3

章 相互作用遮蔽シートを用いたき

裂進展に関する検討

き裂のメカニズムを検討していくため，まず bcc-Fe 完全結晶の中に初期欠陥とし てき裂を導入し，その挙動を検討した．ここでは bcc-Fe 完全結晶内に相互作用を遮 蔽するシートを挿入し，曲率半径 ρ が無限小となる理想き裂を導入して，まずはき 裂進展の様子を観察するとともに局所格子不安定性の観点からき裂進展時の様子を 考察した．

3.1

シミュレーション条件

図 3.1.1 左に模式的に示すように，平板状の周期セル中央に貫通き裂を配置した 系を解析対象とする．き裂前縁方向とき裂面の表記で [001](110)，[11¯2](111) の系 を対象とし，それぞれの結晶方位は x, y, z=[110], [¯110], [001] および x, y, z=[110], [111], [11¯2] である．表 3.1.1 にセル長さ Lx, Ly, Lzおよび総原子数を示した．き裂上 下の原子は互いに力を及ぼさないように制御することでシート状の理想き裂を導入 した．き裂幅 2c は 0.3Lx(x 方向セル長さ) とした．熱揺動による影響を排除するた め，温度は極低温 (0.1K) とし，速度スケーリング法により制御した．10 ps の緩和 計算を行った後，x, z 方向のセル長さを固定して y 方向に引張ひずみを毎ステップ ∆εyy = 1.0× 10−7増加させて引張シミュレーションを行った．数値積分の時間ス テップは 1 fs であるのでひずみ速度に換算すると ˙εyy = 1.0× 108/s である．初期緩 和を除いて，x，z 方向の応力制御は行っていない．

Table 3.1.1 Cell size and number of atoms. model (110) ［111］ Lx[nm] 24.3 24.3 Ly[nm] 32.4 29.8 Lz[nm] 2.0 2.1 Atoms 127,337 122,157 L y Lx Lz 2c=0.3Lx

x

y

z

3.2

シミュレーション結果および考察

3.2.1

応力ひずみ曲線

図 3.2.1 に [001](110) モデルの応力ひずみ曲線を示す．応力ひずみ曲線は引張り開 始から線形に上昇し，εyy = 0.0415 で応力ピークに達して急激に応力が低下し破断 した．

S

tr

es

s,

σ

yy

[

G

P

a]

Applied strain,

εyy

1st peak

ε

yy=0.0415

0 0.02 0.04 0.06

5 10

図 3.2.2 に [11¯2](111) モデルの応力ひずみ線図を示す．引張り開始とともに引張応 力は増加していき εyy = 0.0403 で応力ピークとなった．その後は (110) き裂の応力 ひずみ曲線とは異なり，応力は急激はせずに徐々に減少していく形となった．

S

tr

es

s,

σ

yy

[

G

P

a]

Applied strain,

εyy

1st peak

ε

yy=0.0403

0 0.02 0.04 0.06

5 10

3.2.2

原子応力の分布

図 3.2.3 から図 3.2.6 は，それぞれ [001](110) と [11¯2](111) の応力ピーク時におけ る引張方向の原子応力 σα

yyに関する図である．図 3.2.3 と図 3.2.5 は (a) が各原子を

σα_yyで色付けしたスナップショットであり，色付けの閾値は-10GPa∼80GPa として いる．(b) は格子欠陥を可視化するため可視化ソフト AtomEye の Central Symmetry Parameter(CSP) にて欠陥と判定された原子を示している．(110) き裂ではき裂の上 下斜め方向に高い σα yyが分布している．また (b) を見ると高い σyyα をもつ領域は bcc 構造とは異なった格子欠陥であることが確認できる．(111) き裂は，き裂先端の応 力は (110) のそれよりも低い．CSP による欠陥も殆どはき裂表面近傍にのみ存在す る．図 3.2.4 と図 3.2.6 はき裂面上 (0.5Ly 地点) を y = 0 として，y = 0 から 2d(d は FS ポテンシャルのカットオフ値) の範囲内に存在する原子の σα yy をプロットし た図であり，(a) がき裂先端にとった座標を，(b) が σα yyをプロットしたグラフを示 している．なお，図 (a) の赤い原子は後述する ηα(1) _{< 0 の原子である．(b) は横軸} が原子の x 方向の座標であり，縦軸が原子応力 σα yyの値を示している．同図の赤線 は線形破壊力学による無限体中のき裂の応力分布に，周期条件によって隣接するき 裂の応力場を重ね合わせた式(11)_σ yy(r) = KIC [ 1/√2πr + 1/ √ 2π(d− r)]である．r はき裂先端からの x 方向の距離を示しており，d はき裂先端から隣接周期セルのき 裂先端までの x 方向の距離を示している．破壊じん性値 KIC は等方弾性体を仮定 して KIC = √ 2Eγs/(1− ν2) から求めた．ヤング率 E とポアソン比 ν はそれぞれ E = (C11− C12)(C11+ 2C12)/(C11+ C12)，ν = C12/(C11+ C12) より算出した．(110) と (111) の表面エネルギー γsおよび得られた破壊じん性値 KIC を表 3.2.1 に示す． (111) き裂では応力ピーク点での応力分布は破壊力学による予測によく一致する．一 方，(110) き裂では線形破壊力学の式よりも高い応力を示しており，Griffith の理論 よりも高い応力になってからき裂進展した事がわかる．(110) き裂では応力ピーク前 にき裂から欠陥が生じており，き裂先端の σα yyの分布も大きく変わっていることが この差の一因と考える．

(a) σα

yy snapshot (b) CSP view

Fig. 3.2.3 Snapshots of atoms and defects at the stress-strain peak ([001](110)crack).

0.5L

y

yy

r

Crack tip

2d S tr es s,

σ

yy [ G P a] Position of atoms, r [Å] [001](110) 150 200 250 0 20 40 60 80 (a) model of σα yy distribution (b) σyyα distribution

Fig. 3.2.4 Distribution of atomic stress σ_yyα on y = 0 crack plane at the stress-strain peak ([001](110)crack).

(a) σα_yy snapshot (b) CSP view

Fig. 3.2.5 Snapshots of atoms and defects at the stress-strain peak ([11¯2](111)crack).

0.5L

y

yy

r

Crack tip

2d S tr es s,

σ

yy [ G P a] Position of atoms, r [Å] [112](111) 150 200 250 0 20 40 60 80 (a) model of σα yy distribution (b) σyyα distribution

Fig. 3.2.6 Distribution of atomic stress σα

yy on y = 0 crack plane at the

Table 3.2.1 Young’s modulus, Poisson ratio, surface energies and calculated frac-ture toughness. model [001](110) [11¯2](111) E [GPa] 143.039 ν 0.362 γs [J/m 2 ] 1.703 2.139 KIC [MPa √ m] 0.749 0.839

3.2.3

不安定解析と

AES

の固有値に関する検討

図 3.2.7 は (110) き裂における最小固有値 ηα(1) min の変化であり，破線は ηα(1) = 0 の境 界である．ここで，ηα(1) min の値は各ひずみにおいて全原子の中で最も小さかったもので あり，ある原子の ηα(1)_{をモニターし続けたものではないことに注意されたい．図中に} 示した (a)∼(j) の時点におけるスナップショットを ηα(1) _{< 0 の原子を赤く着色して図} 3.2.8 に示した．[001](110) モデルでは引張り初期は系内に ηα(1) _{< 0 となる原子は存} 在していない．引張り開始からすぐに最小固有値の値は低下して負の原子が現れる． 図 3.2.8(a) に示すように (a)εyy = 0.0031 で相互作用を遮蔽したセル中央部に ηα(1) < 0 となる原子が発生し，引張りに伴って (b)εyy = 0.0033,(c)εyy = 0.0055 の図のように 相互作用を遮蔽した上下部で不安定原子が出現と消失を繰り返す．(d)εyy = 0.0074 以降，き裂先端にのみ ηα(1) _{< 0 の原子が存在するように変化し，(e)ε} yy = 0.0167 で η_minα(1)が最大の負値を示した．その後，(e) → (f) で ηα(1) _{< 0 の原子は消失するが，} (f)εyy = 0.0197 で再びき裂先端に ηα(1) < 0 の原子が現れる．(g)εyy = 0.0219 になる とき裂先端から円弧を描くように結晶内部に ηα(1) _{< 0 の原子が発生し，欠陥の発生} によってバタフライ状に ηα(1) _{< 0 の原子が分布している (図 3.2.8(i))．(110) き裂で} は (i)εyy = 0.0415 の応力ピーク点の後，き裂が (110) 面を不安定に進展して破断する (図 3.2.8(j))．進展時の最小固有値の値はほぼ一定 (図 3.2.7，-200< ηα(1) <-150GPa) で，バタフライ形状も保たれている．完全に破断した後には ηα(1)_{の値は正値に回復} する ((110) 表面原子の値)．

Applied strain,

ε

_yy

M

in

im

u

m

o

f

ei

g

en

v

al

u

e,

η

a (1 )

,

[G

P

a]

η

α(1) (a) (f) (e) (d) (c) (b) (h) (g)

(i) stress peak (j)

0

0.01

0.02

0.03

0.04

0.05

0.06

-400

-300

-200

-100

0

100

(a) εyy = 0.0031 (b) εyy = 0.0033

(c) εyy = 0.0055 (d) εyy = 0.0074

(e) εyy = 0.0167 (f) εyy = 0.0197

(g) εyy = 0.0219 (h) εyy = 0.0267

(i) εyy = 0.0415 (j) εyy = 0.0465

Fig. 3.2.8 Snapshots of atoms during tension in the [001](110)crack. Red circles indicate ηα(1) < 0 atoms.

図 3.2.9 は (111) き裂における最小固有値の変化であり，図 3.2.10 はセル全体の変 形の様子（ただし (a)∼(i) は図 3.2.9 とは対応していない）である．引張り開始から き裂が開口するが，ηα(1) _{< 0 の原子が現れたのは ε} yy = 0.0206 であった．その後， 応力ピークの点の図 3.2.10(c) でも ηα(1) _{< 0 の原子は少なく，(c) → (d) でき裂右端} は水平にわずかに進展し，左端も左下方向へわずかに進展する．(110) き裂とは異な り，き裂先端から転位が射出され，周期境界の下で交差した (図 3.2.10(e) および (f) の黒丸で囲った原子)．その後き裂先端が鈍化し (図 (g))，き裂の上下には先の (110) き裂と同様のバタフライ状の ηα(1) _{< 0 の原子が現れている．その後，き裂右端がジ} グザグに進展して，図 3.2.10(i) で x 周期境界まで達した．この間，図 3.2.9 の ηα(1)_min は約-150GPa で一定である．

応力ピーク前後のき裂部分を拡大して表示したものが図 3.2.11 である．図の (a)∼ (f) は，η_minα(1)の図 3.2.9 に示した点に対応する．図では central symmetry parameter による欠陥もあわせて示している．(110) と同じように引張り初期は不安定原子は 系内に存在せず，εyy = 0.0206 でき裂両端に不安定原子が分布する．図 3.2.9 の応力 ピーク点で ηα(1) min が急減しているが，図 3.2.11(b) の CSP を見るとき裂先端の角部に 欠陥が成長しはじめていることがわかる．図 (b) の ηα(1) _{< 0 の分布を詳細に見ると，} 左の先端には ηα(1) _{< 0 の原子が存在し，右はき裂面に η}α(1) _{< 0 の原子が見られな} い．CSP の形状に違いは見られないが，次の図 (c) では左は斜め下に開口するとと もに積層欠陥が上に発生している．その後き裂先端では積層欠陥が拡大し，右端か らはまず右下方向に転位が射出される ((e)εyy = 0.0412)．さらに右端からは図 (f) の ように下方にも転位が射出された．これらは (112) および (001) のすべり面である． 図 3.2.12 はき裂左端を拡大したものである．(112) 面に開口した際に先端から積層 欠陥を生じている．右端と左端でこのような違いを生じたのは力の遮蔽シートが原 子配置に対して完全に対称となっていなかったためと考える．

Applied strain,

ε

_yy

M

in

im

u

m

o

f

ei

g

en

v

al

u

e,

η

a (1 )

,

[G

P

a]

η

α(1)

_ε

yy=0.0403 stress peak (a) (b) (c) (d) (e) (f)

0

0.02

0.04

0.06

-400

-300

-200

-100

0

100

(a) εyy = 0.0000 (b) εyy = 0.0206 (c) εyy = 0.0403

(d) εyy = 0.0406 (e) εyy = 0.0415 (f) εyy = 0.0420

(g) εyy = 0.0555 (h) εyy = 0.0621 (i) εyy = 0.0669

Fig. 3.2.10 Deformation process during tension in [11¯2](111)crack. Red circles indicate ηα(1) _{< 0 atoms.}

(a) εyy = 0.0206 (b) εyy = 0.0403

(c) εyy = 0.0405 (d) εyy = 0.0406

(a) εyy = 0.0403

(b) εyy = 0.0404

(c) εyy = 0.0405

(d) εyy = 0.0413

Fig. 3.2.12 Snapshots of crack propagation and blunting during tension after stress peak in (111). Left: Red circles indicate ηα(1) < 0 atoms. Right: CSP view.

第

4

章 ねじり粒界モデルに関する検討

前章では，予め欠陥としてき裂を導入したモデルに対してシミュレーションを行っ たが，き裂が発生する状態を議論するためには，き裂がない状態のモデルを検討する 必要がある．今回はき裂の発生源としてねじり粒界を想定し，3 つの CSL(Consident Site Lattice) ねじり粒界モデルを用いて粒界き裂の発生について議論した．

4.1

シミュレーション条件

図 4.1.1 にシミュレーションモデルを示す．図 4.1.1(a) に示すように，セルの上半分 を z 軸 [001] 方向を軸として回転させて x，y 方向で周期性が満たされるサイズで切り出 し，セル中央面にねじり粒界を作成した．回転角度は θ = 36.9◦, θ = 22.6◦, θ = 28.1◦ の 3 つで，それぞれ Σ5, Σ13, Σ17 ねじり粒界となる．その後，図 4.1.1(b) に示すよ うに上下周期境界上の粒界がセル内に入るように周期境界の下で 0.25Lzだけ z 方向 に原子をシフトさせた．表 4.1.1 に系の寸法および総原子数を示す．全方向周期境界 の下で x, y, z 方向の各応力が 0 となるように応力制御しながら 5, 000 [fs] の緩和計算 を行った後，z 方向に引張ひずみを毎ステップ ∆εzz = 1.0× 10−7増加させて引張シ ミュレーションを行った．数値積分の時間ステップは 1 [fs] であるのでひずみ速度に 換算すると ˙εzz = 1.0× 108 [sec−1] である．引張り中の x, y 方向のセル辺長は固定し ている．熱揺動による影響を排除するため，温度は極低温 (0.1[K]) とし速度スケー リング法を用いて制御した．

L z Lx 1/4 L z 0 L y (a) (b) z [001] x [100]

Twin grain boundary Twist

Fig. 4.1.1 Simulation model for twist grain boundaries.

Table 4.1.1 Size, number of atoms and rotation angle. model Σ5 Σ13 Σ17 Lx[nm] 8.6 11.2 9.7 Ly[nm] 8.6 11.2 9.7 Lz[nm] 25.8 33.8 29.2 Atoms 162,000 358,956 235,824 θ[deg] 36.8667 22.6167 28.0667

4.2

シミュレーション結果および考察

4.2.1

Σ5

ねじり粒界の結果

図 4.2.1 に Σ5 モデルの応力ひずみ曲線，第 1 固有値の最小値 ηα(1) min の変化と，緩和 終了後および主要な点におけるスナップショットを示す．スナップショット中の赤 い原子は ηα(1) _{< 0 の原子である．応力ひずみ曲線は ε} zz = 0.0973 と εzz = 0.1415 で ピークを示し，εzz = 0.1470 でほぼ破断した．η α(1) min の変化では，応力ひずみ曲線では 大きな変化がない εzz < 0.0973 の領域で著しい変化を示している．緩和計算終了後 の εzz = 0.0 において，右図に見られるように粒界部分に ηα(1) < 0 の原子が存在する ため，引張り前から負値 (−200GPa 程度) を示している．第一ピーク前の ηα(1) _{< 0} の変化については後述するが，右図のように最初のピークでは粒界近傍に ηα(1) _{< 0} の原子が多く現れ，第二ピークの図では粒内に縞状に ηα(1) _{< 0 の原子が分布してい} る．これは (112) の双晶境界で後で詳細に述べる．破断は粒界部分から生じた．   ηα(1) min の変化の図に示した (a)∼(k) の点および破断後の (l)（範囲外）に対応する ηα(1) < 0 の原子の分布を図 4.2.2 に示す．図では ηα(1) < 0 原子の分布を議論するた めに斜めから見た図や，ηα(1) _{> 0 の原子も表示させた図などを混在させている．図} のセル境界の赤，緑，紫色の直線はそれぞれ x, y, z 軸に対応する．緩和終了後の (a) では上下の粒界部分に ηα(1) _{< 0 となっている不安定原子が層状に存在している．図} 4.2.1 の ηα(1)_min の変化において，まず (a) → (b) は ηα(1)_{の値がほぼ一定で振動してい} る．この時，図 (b) のように ηα(1) _{< 0 の原子層が上下に広がる．(b) → (c) の変化で} は粒界中央部の原子の ηα(1)_{が正に回復する．この時，粒界部は乱れているが bcc 構} 造のままであった．これは CSL 粒界から数層分の原子が引張りを駆動力として原子 の再配置が行われたためと考える．粒界緩和が完了した点が (d) で，ηα(1) min が上昇し ている．(e) は ηα(1) _{> 0 の原子も表示したものである．(d) → (f) では一度回復した} 粒界部分に ηα(1) _{< 0 の原子を生じ，再び (f) → (g) で粒界面の原子の η}α(1)_{は回復し，}

点で中央と下部結晶の内部に ηα(1) _{< 0 の原子が現れ，図 4.2.1 の η}α(1) min の変化でも不 安定に大きく低下している．その後第一応力ピークで粒内に多数の変形双晶が発生 し，応力減少を生じる．図 4.2.3 は初期配置で{112} 面を抽出したものと，第 2 応力 ピーク時の yz 平面のスナップショットを示したものである．粒内の縞は 90◦_回転さ せると一致する．調べた結果，(121) 面であるので第一ピーク後の欠陥は双晶変形に よるものである．  図 4.2.4 はΣ 5 における第 2 応力ピーク後から破断までのスナップショットである． 上の全体図の，中央部の結晶粒と上の粒界近傍を見るとすべり面に沿った結晶内の ηα(1) _{< 0 の縞が別のすべり面のそれと合体している．この後，下に拡大して示した} ように εzz = 0.1432 から上部粒界の ηα(1) < 0 のすべり面と粒界が交差する部分から き裂が発生した．その後は，y の正方向と奥行き方向にき裂が進展し，εzz = 0.1455 のように破断する．

ε

zz

=0.0

ε

zz

=0.1415

2nd stress peak

ε

zz

=0.0973

1st stress peak

ε

zz

=0.1500

Fig. 4.2.1 Stress-strain curve, change in the minimum 1st eigen value η_minα(1) and snapshots of fracture process in Σ5 twist grain boundaries.

(a) εzz = 0.0 (b) εzz = 0.0270 (c) εzz = 0.0383 (d) εzz = 0.0411

(e) εzz = 0.0411 (f) εzz = 0.0668 (g) εzz = 0.0682 (h) εzz = 0.0815

(i) εzz = 0.0960 (j) εzz = 0.1043 (k) εzz = 0.1415 (l) εzz = 0.1470

Fig. 4.2.2 Snapshots of fracture process in Σ5 twist grain boundaries. Red circles indicate ηα(1) < 0 atoms.

θ=79.695

θ

+90[deg]

Twin plane(121)

(a)

(b)

T

win plane(1

12)

y

z

Fig. 4.2.3 Snapshots of slip plane in Σ5 twist grain boundaries. (a) (112) slip plane in initial structure, (b) ηα(1) _{< 0 distribution after the 1st stress}

z

y

ε_{zz =0.1431} ε_{zz =0.1432} ε_{zz =0.1433} ε_{zz =0.1434} ε zz =0.1435 εzz =0.1455 ε zz =0.1415 2nd peak ε_{zz =0.1431}

4.2.2

Σ13

ねじり粒界の結果

図 4.2.5 に Σ13 モデルの応力ひずみ曲線，最小固有値 ηα(1) min の変化とスナップショッ トを示す．応力ひずみ曲線は Σ5 と同じく 2 つのピーク点を持っているが，最小固 有値 ηα(1) min の変化は Σ5 と比較して著しい変化はなく，振動しながら低下している． スナップショットでは粒界部分に ηα(1) _{< 0 原子が分布 (ε} zz = 0.0) し，応力ひずみ 曲線の第 1 ピークで広がっているのは Σ5 と同じである．第一ピーク後に粒内に縞 状の ηα(1) _{< 0 の分布を生じるのも Σ5 と同じであるが，Σ5 と異なり Σ13 では縞は} 上下の粒内で平行となっている．第二ピーク後，εzz = 0.1530 で粒界から破断した． 図 4.2.6 は図 4.2.5 の ηα(1) min の (a)∼(l) に対応させて ηα(1) < 0 の分布を見たものであ る．Σ5 と異なり，粒界部分の回復などはなく粒界近傍の ηα(1) _{< 0 の層が増加する．} 応力ピーク直前の ηα(1) _{< 0 の分布 (e) と (f) の CSP 表示を見ると欠陥だけでなく上} 下の原子層も ηα(1) _{< 0 となっている．応力ピーク後の (g) → (h) では下部粒界から} ηα(1) _{< 0 の集団が面状に成長している．Σ5 モデルでは，結晶内に η}α(1) _{< 0 の原子が} 数個現れた後に粒内に一気に双晶が形成されたが，Σ13 では転位ループのように粒 界から伝ぱしている．しかしながら，応力ひずみ曲線の谷に相当する εzz = 0.09994 の (i),(j) 図では，縞状に不安定原子が分布している (j) に対して (k) の CSP で多くの 原子が欠陥原子と判定されており，1 つの転位がすべり面を横断したような変形で はない．  図 4.2.7 は応力第 2 ピークにおける yz 平面のスナップショットを初期配置におけ る_{{112} すべり面と比較して示したものである．η}α(1) _{< 0 の縞の向きは Σ5 と異な} り，_{{110} すべり面に一致する．上下で模様の差が少ないのはねじり角 θ が小さい} ためと考えられる．  図 4.2.8 にΣ 13 におけるき裂進展の様子を示す．左の全体スナップショットを見 ると第 2 応力ピーク時 (εzz = 0.1435) に中央結晶内に発生した縞模様が εzz = 0.1495

復範囲が大きいことと，先の図 4.2.6(k) の CSP から面間隔が大きくひずんだ部分が ηα(1) _{< 0 の縞状の構造と推測する．下の拡大図を見ると斜めの η}α(1) _{< 0 の縞によっ}

て粒界がジグザグに波打った状態となっている．εzz = 0.1496 で上部粒界の周期境

ε

zz

=0.0

ε

zz

=0.1435

2nd stress peak

ε

zz

=0.0909

1st stress peak

ε

zz

=0.1530

Fig. 4.2.5 Stress-strain curve, change in the minimum 1st eigen value η_minα(1) and snapshots of fracture process in Σ13 twist grain boundaries.

(a) εzz = 0.0 (b) εzz = 0.0316 (c) εzz = 0.0525 (d) εzz = 0.0690

(e) εzz = 0.0909 (f) εzz = 0.909(CSP) (g) εzz = 0.0910 (h) εzz = 0.0913

(i) εzz = 0.0925 (j) εzz = 0.0994 (k) εzz = 0.994(CSP) (l) εzz = 0.1530

θ = 80.049

θ

θ'

θ' = 69.964

slip plane{110}

(a)

(b)

twin plane{1

12}

z

y

Fig. 4.2.7 Snapshots of slip plane in Σ13 twist grain boundaries. (a)slip plane model, (b) simulation model with yz view in 2nd stress peak.

εzz =0.1495

z

y

ε_zz =0.1495 ε_zz =0.1496 ε_zz =0.1497 ε_zz =0.1498 ε_{zz =0.1499} ε_{zz =0.1518} εzz =0.1435 2nd peak

4.2.3

Σ17

ねじり粒界の結果

図 4.2.9 に Σ17 モデルの応力ひずみ曲線および最小固有値 ηα(1) min の変化とスナップ ショットを示す．やはり応力ひずみ曲線は大きな 2 つのピークを示した後破断した が，最小固有値 ηα(1) min は Σ5 と同じく複雑な変化を第 1 ピーク前に示している．何よ り特徴的なのは第 1 応力ピークでは粒界部以外の原子が ηα(1) _{< 0 となっていること} である (右上スナップショット)．このように，結晶全体が ηα(1) _{< 0 となるときは，} 転位のような局所的な変化よりも相変態のような集団的な変形を生じやすいことを 報告している(8)_{．その後，第 1 ピーク直後に η}α(1) min は急減し，粒内に Σ13 と同様の 縞模様を残して回復する．  図 4.2.10 に ηα(1) min の変化の (a)∼(i) に対応させた図を示す．εzz = 0.0 では系内に ηα(1) < 0 の原子は存在しない．(b) → (c)，(d) → (e) などの変化は Σ5 の変化と似てお り，粒界近傍の構造緩和と考えられる．(f) の時点ではまだ ηα(1) _{< 0 の原子は粒界近} 傍のみだが，(g) の第 1 ピークでは図 4.2.9 に示したように粒内すべてが ηα(1) _{< 0 と} なる．この時を CSP による可視化をしたのが (g) で粒内に欠陥原子は認められない． その後の応力急減をもたらす変形によって粒内の多くの原子は回復して ηα(1) _{> 0 と} なるが，図 4.2.9 の ηα(1) min の (f) → (g) → (h) の変化からわかるように縞の部分の ηα(1) はより大きな負値となる． 図 4.2.11 は応力第 2 ピークにおける yz 平面のスナップショットである．同図 (a), (b) の青およびピンク色の線が{112} すべり面を yz 平面に投影したものである．同 図 (b) より，ねじっていない中央結晶の縞は青線に一致せず，{110} すべり面を投影 したに一致する．しかし，ねじった上下の結晶内は_{{112} すべり面を投影したピン} ク色の線に一致しており，上下の粒内で別の変形を生じたことが推測される．  図 4.2.12 は破断した下の粒界を側面から見たものである．Σ17 では図 4.2.9 の応 力ーひずみ曲線のように小さな第 2 ピークの後，第 3 ピークで粒界をはさんで上下

ε

zz

= 0.0

ε

zz

=0.1535

3rd stress peak

ε

zz

=0.1037

1st stress peak

ε

zz

=0.1595

Fig. 4.2.9 Stress-strain curve, change in the minimum 1st eigen value η_minα(1) and snapshots of fracture process in Σ17 twist grain boundaries.

(a) εzz = 0.0 (b) εzz = 0.0400 (c) εzz = 0.0647

(d) εzz = 0.0649 (e) εzz = 0.0870 (f) εzz = 0.0965

θ = 78.9174

θ

θ'

θ' = 71.8351

slip plane{110}

slip plane{112}

θ+90[deg]

slip plane

{1

12}

slip

plane{1

12}

x

z

Fig. 4.2.11 Snapshots of slip plane in Σ17 twist grain boundaries. (a)slip plane model, (b) simulation model with yz view in 2nd stress peak.

z

y

ε zz =0.1539 εzz =0.1540 ε_{zz =0.1541} ε_{zz =0.1542} ε zz =0.1543 εzz =0.1575 εzz =0.1540 3rd stress peak

第

5

章 円孔の応力集中部からのき裂発

生に関する検討

前章では粒界き裂の発生を検討したが，どのモデルもき裂発生前にまず粒界から 結晶内部への塑性変形を生じ，粒内に縞状の不安定原子が分布した．粒界近傍の不 安定原子の割合などを調べたが，き裂発生箇所を特徴づけるような有意な結果を見 出すことは難しかった．そこで今回は奥行き方向のセル長さを小さくすることで平 面的な固有ベクトルの観察を実現するとともに，き裂の発生源として応力集中部（円 孔）を導入し，応力集中部からのき裂発生について検討した．

5.1

シミュレーション条件

図 5.1.1 に示すような bcc-Fe の完全結晶薄板セルの中央に，0.3Lxの大きさとなる 円孔を設けることで応力集中部を形成した．x, y, z を bcc 構造の [110]，[¯110]，[001] と する [110] モデルと，[110]，[111]，[11¯2] とする [111] モデルの 2 つを対象とした．系の 寸法および原子数を表 5.1.1 に示す．全方向周期境界条件下で x, y, z 方向の各応力が 0 となるように応力制御しながら 10, 000 [fs] の緩和計算を行った後，x, z 方向のセル辺 長は固定する条件下で y 方向に引張ひずみを毎ステップ ∆εyy = 1.0×10−7増加させて 引張シミュレーションを行った．数値積分の時間ステップは 1 [fs] = 1.0× 10−15_[sec] であるのでひずみ速度に換算すると ˙εyy = 1.0× 108 [sec−1] である．熱揺動による 影響を排除するため温度は極低温 (0.1[K]) とし，速度スケーリング法を用いて制御 した．

Table 5.1.1 Cell size and number of atoms. model ［110］ ［111］ Lx[nm] 24.3 24.3 Ly[nm] 32.4 29.8 Lz[nm] 2.0 2.1 Atoms 127,337 122,157

L

y

Lz

Lx

y x

5.2

シミュレーション結果および考察

5.2.1

応力ひずみ線図および最小固有値の変化

図 5.2.1 および図 5.2.2 に応力ひずみ曲線，および，全ての原子の中で最小の ηα(1) min の値をプロットしたものを示す．左目盛りに応力，右目盛りに ηα(1) min の値をとってい る．水平破線が固有値 0 の境界を示している．応力はいずれもほぼ線形に上昇した 後，ピークを示し応力急減しているが，詳細に見ると [110] モデルでは矢印で示した εyy = 0.0375 で小さな応力減少を示しており，その後の応力ひずみ曲線にも小さな振 動が認められる．εyy = 0.0375 の直前に AES の最小固有値 η α(1) min の値が急減し，矢印 で示した点以降で大きく揺らいでいる．εyy = 0.0562 の応力ピーク点では，後述する ようにき裂が発生・進展した． [111] モデルの場合，最小固有値が負になるのは応力

S

tr

es

s,

σ

yy

G

P

a

strain,

εyy

M

in

im

u

m

o

f e

ig

en

v

al

u

e,

η

a (1 )

, G

P

a

η

α(1)

σ

yy [110] model 0 0.02 0.04 0.06 5 10 15 -200 -100 0

Fig. 5.2.1 Stress-strain curve and change in the minimum 1st eigen value ηα(1) under [110] tension.

S

tr

es

s,

σ

yy

G

P

a

strain,

εyy

M

in

im

u

m

o

f e

ig

en

v

al

u

e,

η

a (1 )

, G

P

a

η

α(1)

σ

yy [111] model 0 0.02 0.04 0.06 5 10 15 -400 -300 -200 -100 0

Fig. 5.2.2 Stress-strain curve and change in the minimum 1st eigen value ηα(1)

ピークの直前であり，その前の応力ひずみ曲線には [110] モデルのような変化は認め にくい．応力ピーク点では転位射出が複数回生じ，それに対応して応力が段階的に 減少しており，き裂は発生していない．図 5.2.3 および図 5.2.4 は [110]，[111] の応力 ピーク前後のスナップショットである．AES の第 1 固有値 ηα(1)_{が負の原子を赤く着}

色している．[110] モデルの図 5.2.3(a) において，εyy = 0.0562 で円孔を中心に斜方に

不安定原子が分布している部分は AtomEye の Central Symmetry Parameter(CSP) でも積層欠陥と判定されており，先の応力ひずみ曲線で変化が現れた時に欠陥がバ タフライ状に成長したことを確認している．ピーク点以降，円孔赤道部からき裂が発 生し図 5.2.3(b) の図のように水平方向へ不安定に進展して破断した．図 5.2.4 の [111] モデルでは εyy = 0.0511 の応力ピーク時に [110] モデルと同じく円孔斜方に ηα(1) < 0 の原子がわずかに存在している．しかし，[110] と異なり円孔赤道部には ηα(1) _{< 0 の} 原子は発生していない．εyy = 0.0514 近傍で円孔斜方の不安定原子群の根本から赤 道部へ向かって不安定原子（転位）が射出された．(111) 面での角度からこの転位は (112) 面と考えられる．図 (c) では ηα(1) < 0 の原子が交差し 1 つに見えるが，CSP に より欠陥を可視化すると 2 つの転位芯が確認できる．その後 εyy = 0.0527 で x 方向 周期境界に達して反対側の円孔表面から射出された転位と交差し停滞した．このよ うな切り欠きからの転位射出や他の転位との交差など複雑な応答を繰り返して塑性 変形し，εyy = 0.0800∼0.0900 の変形後期に図 (e) のような小さなき裂を生じるもの の，[110] と異なりき裂の不安定成長は生じずに鈍化と進展を繰り返した．

cracking

(a) εyy = 0.0562 (b) εyy = 0.0570

Fig. 5.2.3 Snapshots before and after the stress-strain peak ([110] tension).

con

emission of dislocation

cracking

CSP view

5.2.2

固有ベクトルによる検討

ηα(1) < 0 の原子について，対応する固有ベクトル{∆εi} の成分 ∆εxx, ∆εyy, ∆εzz, ∆εyz, ∆εzx, ∆εxy からひずみテンソル [∆εij] を求め，その主軸である固有ベクトル x を求めた．固有値 ε1 < ε2 < ε3に対応する固有ベクトル x1, x2, x3のうち x3につ いて，ε3の大きさに応じてスケーリングし xy 平面に投影したものを図 5.2.5 および 図 5.2.6 に示す．詳細観察のため原子座標の出力間隔を細かくしてとりなおしたの で図 5.2.3 および図 5.2.4 で示したスナップショットのひずみとはわずかに異なる． 図 5.2.5(a) は応力ひずみ線図のピーク時で，き裂開口する直前である．き裂が発生 する図 (b) では左のグラフ上に (010) 面を破線で示した．き裂開口のする先端では ηα(1) < 0 の原子の変形モードは (010) 面に対して垂直であることが分かる．き裂進 展中の図 (c) でも不安定モードは (010) 面に垂直または平行であり，(110) 面で開口 しているが開口した側の結晶方位は (010) 面に沿う方向に傾いている．図 (c) にこの 変形モードの方向を矢印で示した． 図 5.2.6(a) は [111] 引張りで，第 1 応力ピーク後に ηα(1)_{がパルス状に大きな負値} を示したときの図である．同図左には_{{112} すべり面を破線で示した．(a) → (b) に} かけて円孔斜方から不安定原子の射出が見られる．射出された不安定原子の x3には 規則性は見られないが，転位射出元である円孔斜方の不安定原子は (112) すべり面 に対して，垂直もしくは水平な x3のどちらかを持っている（青丸印）．図 (c) では 円孔斜方から不安定原子が射出され転位となるが，射出された転位は図 (b) と同じ く_{{112} すべり面に対して垂直か水平なベクトルが存在している．転位はすべり面} 上を動いており，不安定原子は転位によって結晶格子の形状を乱されるため，垂直 方向のベクトルはパイエルスポテンシャルに相当する原子面間方向のモードと考え られる．

(a) εyy = 0.05620 [Å] [Å]

x

₃ (b) εyy = 0.05624 [Å] [Å]

x

₃ (c) εyy = 0.05670

(a) εyy = 0.05120 [A] [A] x₃ {112} slip plane (b) εyy = 0.05130 [A] [A] x₃ {112} slip plane (c) εyy = 0.05140

Fig. 5.2.6 Principal axis x3of unstable deformation mode of [∆εij] under the [111]

第

6

章 結言

原子弾性剛性係数 Bα ij ≡ ∆σαi/∆εjの正値性によって変形・破壊を統一的に議論す る研究の一環として，bcc-Fe のき裂を対象とした様々な分子動力学シミュレーショ ンを行い，Bα ij の固有値 ηαεiの変化や分布，固有ベクトルなどを議論した．   2 章では，分子動力学法および EAM ポテンシャルについて概説し，弾性剛性係 数と格子不安定性について説明した．また EAM における原子応力 σα i や原子弾性剛 性係数 Bα ijの導出とそれから得られる Bijα について示した．   3 章では，原子間相互作用を遮蔽するシートによるき裂の引張りシミュレーションを 行った．[001](110) き裂では応力ピーク後にき裂が水平に進展し破断した．[11¯2](111) き裂では，わずかに進展したものの，その直後に転位射出して鈍化した．き裂面の 原子応力と線形破壊力学による応力分布を比較した結果，応力ピーク前に欠陥を生 じなかった [11¯2](111) き裂は線形破壊力学の式に一致したが，[001](110) き裂では応 力ピークより前にき裂から欠陥原子が発生し，破壊力学による予測よりも原子応力 の方が高い値をとった．全ての原子の AES の第 1 固有値の ηα(1)_{のうち最小のもの} (η_minα(1)) の変化を調べたところ，どちらのき裂でも原子構造の組み替え，すなわち欠 陥を発生する直前に ηα(1) min の最小値が急減し，-200GPa< η α(1) min <-150GPa で変形が進 行した．   4 章ではき裂の発生を議論するため，き裂発生箇所としてねじり粒界を想定し， Σ5, Σ13, Σ17 の 3 つの CSL(Consident Site Lattice) ねじり粒界を周期境界によって ラミネート状に有するバルクの引張りシミュレーションを行った．いずれも応力ひ