原子弾性剛性係数に基づく局所不安定性解析：fccニッケル中のき裂先端の力学特性評価

修士論文

原子弾性剛性係数に基づく局所不安定性解析

：

fcc

ニッケル中のき裂先端の力学特性評価

指導教員

:

田中克志

 

平成

28

年

2

月

神戸大学工学研究科

博士課程前期課程 機械工学専攻

147T327T

片山 寛

Master thesis

Local lattice instability analysis

based on atomic elastic stiffness:

study on crack tip in fcc nickel.

Hiroshi KATAYAMA

February 2016

Department of Mechanical Engineering,

Graduate School of Engineering,

要約

本論文では，モードI型き裂を有するfcc-Niを対象とし，分子動力学法を用いて引張シ ミュレーションを行い，き裂進展挙動と原子弾性剛性係数B_ijα の行列式detB_ijα，ならび にその固有値ηα の関係について検討を行った．まず[001](010)，[¯110](001)，[¯1¯12](111) の三種類の結晶方位で平板状の周期セル中央に貫通き裂を有する系に対し引張シミュレー ションを行った．いずれのモデルも無負荷平衡状態ではdetB_ijα < 0の原子は存在して おらず，き裂近傍に detBα ij < 0となる原子が生じた後，き裂先端から欠陥が発生し応 力-ひずみ線図に非線形性が現れた．また固有値の変化を調べた結果，構造変化が起こる 点(欠陥が発生する点)ではηα(1)の最小値がほぼ一定の負の値([001](010)，[¯1¯12](111) は-40GPa，[¯110](001)は-20GPa）となっていた．次に周期セルの中央と四隅に貫通き裂 を配置したダブルき裂モデルでの解析を行った．欠陥生成時のdetBα_ij < 0原子の数はい ずれも単一き裂モデルのほぼ二倍であり，き裂先端の不安定領域は一定の寸法である可能 性が示唆された．次に周期セルの厚さ方向を1.5倍とし，厚さ方向自由度を変えたモデル で解析を行った．厚さ方向の自由度を増やしたことによって変形モードがわずかに変化し たが，detB_ijα < 0の原子が初めて現れるひずみにはどのモデルも変化がなく，その原子 数も1.5倍([001](010)のみ1.4倍)となっていることから，き裂先端にdetB_ijα < 0の原 子が発生するときの力学条件は厚さ方向の自由度の影響はない．全ての解析を総合して， 結晶方位による差はあるが，き裂先端に初めてdetB_ijα < 0 が出現するときの局所力学状 態は常に同じ(detBα ij < 0の原子列＝不安定領域が一定)であり，またその後き裂先端か ら転位射出しき裂が鈍化・開口するときのBα ij の第一固有値が一定の値を示すという結果 を得た．

Summary

In order to discuss the local stability at a crack tip, we performed molecular dy-namics (MD) simulations on mode I through crack in fcc nickel and evaluate the determinant of atomic elastic stiffness B_ijα, and its eigenvalue ηα. First, we per-formed tensile simulations on single crack models with three crystal orientation of the [001](010)，[¯110](001)，and [¯1¯12](111), respectively. There is no detB_ijα < 0 atoms

in the no load equilibrium, and detBα

ij < 0 atoms appear at the crack tip at the

strain 0.035，0.042，and 0.026.Then inelastic deformation takes place at the tip and

the stress-strain curves show nonlinearity at the point. The minimum 1st ηα(1) of

B_ijα becomes negative and shows almost constant value,-40GPa for [001](010), and

-20GPa for [¯110](001), during crack propagation. Next, we simulate with the double crack model of center and corner crack in the periodic cell in order to discuss the effect

of different crack array. The number of detBα

ij < 0 atoms is simply doubled at the

defect nucleation, so that we conclude that the size of the unstable region is almost constant, despite of the different force field. Finally, we analyses thicker model by 1.5 times to change the flexibility of the thickness direction. Deformation mode is slightly changed only in the [¯1¯12](111) crack; however, the strain at the first emerging point of detB_ijα < 0 atoms is identical to the precious ones. The number of detB_ijα < 0 atoms is just 1.5 times (the [001](010) model 1.4 times) larger than that of the base models. So the flexibility of thickness direction doesn’t affect the mechianical condition for the emergence of detB_ijα < 0 atoms at the crack tip. As a summary, we observed the local criteria for the emergence of detB_ijα < 0 atoms at the crack tip, the 1st eigenvalue of detB_ijα for the subsequent defect nucleation which leads the crack propagation.

目次

1 緒論 1 2 解析手法の基礎 4 2.1 分子動力学法 . . . 4 2.2 速度スケーリング法 . . . 5 2.3 原子埋め込み法ポテンシャル . . . 7 2.4 高速化手法 . . . 8 2.5 原子弾性剛性係数 . . . 9 3 モードI型き裂の引張シミュレーション 15 3.1 シミュレーション条件 . . . 15 3.2 シミュレーション結果および考察 . . . 17 4 ダブルき裂モデルでのシミュレーション 29 4.1 シミュレーション条件 . . . 29 4.2 シミュレーション結果および考察 . . . 30 5 厚さ方向自由度の影響 38 5.1 シミュレーション条件 . . . 38 5.2 シミュレーション結果および考察 . . . 38 6 結論 46 参考文献 48 謝辞 50

1

緒論

分子動力学シミュレーションは材料の特性を個々の原子の運動にまで遡って解析する手 法であり，材料科学だけでなく，高分子化学や生体化学など様々な分野で活用されてい る．計算機が飛躍的に性能の向上，およびナノレベルでの材料の変形・破壊実験が可能と なった今日，ナノ∼サブミクロン領域の現象を疑似的に観察する手法として分子動力学シ ミュレーションは広く用いられるようになってきた．しかしながら，分子動力学シミュ レーションの中で見出した原子レベルでの変形メカニズムについて，どのような指標を用 いて普遍的・定量的な議論に結び付けるかということは依然として問題である．実在の材 料の挙動をできるだけ精密に再現するポテンシャル関数の開発に膨大な努力がなされてい るが，「似ているがそのものではない」モデル材料の挙動であるため，局所エネルギーや 応力などのしきい値を求めても意味がない．そのため応力拡大係数などの連続体理論を指 標とした議論に頼りがちだが，き裂先端などの原子構造に起因したわずかな違いなどは議 論できず，モデル材料における連続体理論の適用限界を示すのみとなる． 著者らのグループでは，局所変形抵抗を表す物理量として原子弾性剛性係数(Atomic

Elastic Stiffness, AES)に着目し，その正値性を基準とすることで，用いたモデル材料や

ポテンシャルエネルギーによらない普遍的な支配力学法則の解明を目的とした議論を行っ てきた1)．AESによる検討は，fccのNiやbcc構造のFe，アモルファス金属やダイヤモ ンド構造を持つSiなどに広く展開してきた．1)–4) 本研究で対象とするfcc構造のNiにつ いては，単結晶ナノワイヤの引張シミュレーションより，表面から部分転位が発生・移動 する際にはAESが負の原子が転位前縁に現れることを明らかにしている1)．アモルファ ス構造への適用では，AESの正値性によって局所クラスタの崩壊及び再構成が理解でき

ることを報告している2)．Siの検討では，表面や界面による内部不均一性に着目した検討 を行い，これらの欠陥部分のAESは無負荷平衡状態で負または0に近い値を示し，欠陥 から離れたバルクの部分と大きな差があるものの，外力を与えるとバルク部分のAESが それに近づく「系の均一化」を生じていることを報告している3)．また，bcc鉄中のらせ ん転位に関する検討では，転位芯原子のAESの行列式の値は無負荷平衡状態では負とな らず，外力を受けて転位が移動する際に初めて負となることなどを明らかにしている4)． 分子動力学シミュレーションによるき裂の解析は古くから行われてきた．計算機能力が 十分でなかった1970年代に既に，き裂の開口時の原子結合切断を簡略化してモデル化す ることで，格子力学的な視点からき裂の進展を議論した検討がなされている5)．Gumbsch らはき裂の進展を局所格子の破断と捉えたモデルを提案し，第一原理計算を用いた精密な 評価を行っている6)．また，Shimadaらは系の全自由度に対してヘッシアンを計算する という厳密な手法を用いて，Siき裂についてき裂進展の解析を行い，き裂先端の一格子程 度の領域のみが不安定となっていることを明らかにした7)．しかし系のヘッシアンを直接 扱う方法では計算量が膨大となるため，100万原子程度の大きな系では直接解くことがで きず，き裂近傍の一部をとりだして弾性体中の境界条件を新たに与えて解いている．一方 著者らのグループで提案しているAES解析は、後で述べるように6× 6のマトリックス で表されるので，全原子について固有値解析を行うことが容易である．そのため大規模な 系での解析を行うことが可能であり，また熱搖動の影響やき裂の鈍化についても考慮する ことができる．き裂を対象としたAES解析はこれまでにbcc-鉄,シリコンについて進め てきた8), 9)．bcc鉄での解析では，き裂が進展する際に塑性変形がほとんど発生しない結 晶方位([001](110),[¯1¯12](111))では，き裂の大きさや数，配置の仕方に関わらずき裂先端 に現れる負のAES原子数が一定であることを明らかにしている9)．またSiの解析では，

転位の射出を妨げSiのぜい性的な挙動をもたらしている可能性や，き裂進展直前に大き な負の第一固有値を生じる前駆挙動があること，またその固有ベクトルは[001](010)き裂 はモードI，[001](110)き裂はモードII，[112](111)き裂ではモードIIIの変形モードで あることなどを明らかにしている8)． 本研究では，fcc金属のNiをモデル材料として同様のモードI型き裂の引張シミュレー ションを行い，き裂進展/鈍化時のAESを明らかにすることを目的とする．以下に各章の 概略を示す．第2章ではMD法の基礎理論を述べ，用いた原子間ポテンシャルにおける 原子弾性剛性係数の定式化を行う．第3章ではモードI型き裂の引張シミュレーションを 行い，その変形挙動をAES解析を用いて検討を行う．第4章ではき裂を複数導入するこ とでき裂配置の影響について検討を行う．第5章ではシミュレーションセルの厚さ方向の 原子数を変えたモデルを用いて解析を行い，厚さ方向の自由度の影響について検討する． 第6章では得られた結果を総括する．

2

解析手法の基礎

2.1

分子動力学法

分子動力学法(molecular dynamics method，略してMD法)は，系を構成する各粒子

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

であり，式(2.1)と式(2.5)を式(2.3)と式(2.4)に代入すると， rα(t + ∆t) = rα(t) + ∆tvα(t) + (∆t) 2 2 Fα(t) mα + O ( (∆t)3 ) (2.6) rα(t− ∆t) = rα(t)− ∆tvα(t) + (∆t) 2 2 Fα(t) mα + O ( (∆t)3 ) (2.7) となる．両式の和と差をとると， rα(t + ∆t) + rα(t− ∆t) = 2rα(t) + (∆t)2 F α (t) mα + O ( (∆t)4 ) (2.8) rα(t + ∆t)− rα(t− ∆t) = 2∆tvα(t) + O ( (∆t)3 ) (2.9) が得られる．これより，時刻t + ∆tでの位置ベクトルとtでの速度は rα(t + ∆t) = 2rα(t)− rα(t− ∆t) + (∆t)2 F α (t) mα + O ( (∆t)4 ) (2.10) vα(t) = 1 2∆t {r α (t + ∆t)− rα(t− ∆t)} + O ( (∆t)2 ) (2.11) と求められる．t + ∆tでの座標を求めるには2つの時刻tとt− ∆tでの座標が必要であ る．初期の計算(t = 0)では，t = ∆tでの座標rα(∆t)は式(2.6)と初速度から求められ る．ri(∆t)とri(0)が既知であれば，式(2.10)を繰り返し適用することにより各粒子の 座標を求められる．

2.2

速度スケーリング法

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

mα :粒子αの質量 vα_i :温度T での粒子αの速度 kB : Boltzmann 定数 = 1.38× 10−23[J/K] 目標の温度 T0 における原子 α の速度を viα0 とおくと v α i0 は式(2.13) のように表さ れる．  vα_i 0 = ( 3kBT0 mα )0.5 (2.13) 同様に，温度 T の時の原子 α の速度は式(2.14) のように表される．  vα_i = ( 3kBT mα )0.5 (2.14) よって，式(2.13)と式(2.14)より以下の式が得られる． vα i0 vα i = ( T0 T )0.5 (2.15) つまり，系の温度をT から T0 にするには，式(2.15)の右辺を現在の速度に掛けてやれば よい．ただ，これだけでは原子配置に反映されないので，Verlet法における∆rα_i(t + ∆t) （式2.16）を√T0/T ∆rαi(t + ∆t)と置き換える必要がある． ∆rα_i(t + ∆t) = rα_i(t + ∆t)− rα_i(t) = rα_i(t)− rα_i(t− ∆t) + (∆t)2F α i (t) mα (2.16) 平衡状態では，能勢の方法10) など外部との熱のやりとりをする変数を考慮した拡張系の 分子動力学法によって得られるカノニカルアンサンブルに一致することが示されている．

2.3

原子埋め込み法ポテンシャル

分子シミュレーションにおいて，系の振る舞いを精度良く再現するためには分子に作用 するポテンシャルエネルギーが重要となる．量子力学に基づき，電子や原子核のハミルト ニアンから系のポテンシャルエネルギーを精密に求める第一原理計算も試みられている が，計算量が極めて膨大となるため扱える粒子数が限られてしまう．そこで多数の原子の 挙動を求める場合には原子間相互作用を簡略評価する原子間ポテンシャルが通常用いられ る．これまでに数々のポテンシャル関数が提案され用いられてきたが，解析の目的に対し て適切なポテンシャルの設定を誤ると，いかに計算を重ねたところで無意味な検討となっ てしまうため，ポテンシャル関数の選択は慎重に行わなければならない．

原子埋め込み法(Embedded atom method)12)13)はDawとBaskesによって Stotと

Zarembaの擬原子法 14) の考えを用いて提案されたポテンシャル関数で，金属中の多体 効果を精度良く再現することから広く用いられている．EAMでは密度汎関数理論に基づ き，まず金属材料における系全体のポテンシャルエネルギーEtot は原子を価電子雲中に 埋め込むエネルギーと原子間の２体間相互作用の和で与えられるとする．さらに，埋め込 みエネルギーは埋め込む位置の電子密度にのみ依存すると仮定することによって，系全体 のエネルギーは次式のように表される． Etot = 1 2 ∑ i,j(̸=i) Vij(rij) + ∑ i Fi(ρi) (2.17) ここで，ρiは，原子iの位置における電子密度，Fi (ρi)は電子密度ρi の位置に原子を 埋め込むエネルギー，Vij(rij)は距離rij 離れた原子iとj のクーロン相互作用である．

電子密度ρα は周囲の原子βからの寄与ϕ(rαβ)の重ね合わせで与えられると仮定し ρα = neighbor_∑ β(̸=α) ϕ(rαβ) (2.18) で評価する．本解析ではVoterがNi,Alに対し弾性定数，格子定数，原子エネルギー等へ のフィッティングを行ったEAMポテンシャルを使用する.15) ここでの二体間ポテンシャ ルV (rαβ)，二体間関数ϕα，埋め込み関数F (ρα)はすべてスプラインデータで提案され ており，数式としては与えられていない．本研究においてはこれらを3次のスプライン関 数にフィッティングすることにより使用している．

2.4

高速化手法

原子論に立脚したシミュレーションでは原子間相互作用の計算に大部分の負荷が集中す る．しかし実際の結晶中では近接原子の遮蔽効果により，離れた原子からの相互作用はほ とんど及ばない．そのため，相互作用打ち切り(カットオフ)半径rc を導入し，その半径 内の原子からの寄与のみを考慮すればよい．例えば2体ポテンシャルを用いる場合，系 の全粒子数をN とすると系の内部だけでも組み合わせはN (N− 1)通り存在し，粒子数 が増加することで飛躍的に計算コストが増えることがわかる．この計算コストを減らす方 法として領域分割法を用いた．領域分割法では系を間隔がrc 程度である格子状に分割し， 対象とする原子の存在する領域(図2.1着色部)および隣接する領域(図2.1斜線部)の原 子の内，距離がrc 以下であるものを探索する.    

   

x

y

0

bx

by

  Fig.2.1 Schematic of domain decomposition.

     

2.5

原子弾性剛性係数

2.5.1 弾性剛性係数と格子不安定性 応力 σij および弾性係数 Cijkl は，等温過程では σij = 1 V ( ∂F ∂ηij ) , Cijkl = 1 V ( ∂2F ∂ηij∂ηkl ) (2.19) と定義される16)．ここで，F はHelmholtzの自由エネルギー(断熱過程では内部エネル ギーU )，V は結晶の体積，ηij は平衡状態(無負荷とは限らない) からの仮想的な微小ひ ずみである．一方，無負荷平衡状態を基準とするひずみεij と応力σij の関係は，2つの

平衡状態間の変形を考えて導出される次の弾性剛性係数によって表される16)． Bijkl ≡ ( ∂σij ∂εkl ) = Cijkl + (σilδjk + σjlδik + σikδjl+ σjkδil− 2σijδkl)/2 (2.20) ここで，δij はクロネッカーのデルタである．Bijkl は非線形弾性における応力-ひず み関係の勾配を表すことから，Wang, Yip らは，ひずみの対称性を考慮したテンソル

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

る．Bijkl の最小固有値が負であれば，対応するひずみに対して応力が負になる自発的変

2.5.2 原子系への拡張 局所の安定性を評価するための原子弾性剛性係数B_ijklα の算出に必要な原子応力σ_ijα な らびに原子弾性係数C_ijklα は，各原子周りの微小ひずみに対するポテンシャルエネルギー の1次，2次変化量として導出される． 2.5.3 原子応力 簡単のため，結晶の内部エネルギー U が系全体のポテンシャルエネルギーEtot に等し いとする．このとき，応力は平衡状態からの微小ひずみ η に対するポテンシャルエネル ギーの単位体積当たりの変化として与えられる16)． σij = 1 V ∂Etot ∂ηij (2.21) ここで，V は平衡状態における系の体積であり，下付添字のローマ文字はテンソルのデカ ルト座標成分を表す．(2.21)式の微分を求めるため，平衡状態からの仮想的な均一変形を 考える．結晶内の原子 α の位置ベクトルは仮想変形のヤコビ行列 J によって rα = J ¯rα (2.22) と変化する．ここで，「¯」は仮想ひずみによる変形前の値を示す．これより，原子 α と 原子 β の間の距離 rαβ には (rαβ)2 = ¯r_iαβGij¯rαβj (2.23)

なる関係が成立する．ただし，Gij = JkiJkj である．仮想変形のLagrangeひずみテン ソル ηij は ηij = 1 2 [ Gij − δij ] (2.24) であり，その微小量 dηij = 1 2dGij (2.25) と式(2.23)の関係から次の関係が得られる． ∂rαβ ∂ηij = r¯ αβ i ¯r αβ j rαβ (2.26) これよりEAMポテンシャルEtot における応力は次式で評価される σij = 1 V ∂Etot ∂ηij = 1 V ( 1 2 N ∑ α N ∑ β(̸=α) ∂rαβ ∂ηij ∂Etot ∂rαβ ) = 1 V N ∑ α N ∑ β(̸=α) ( F′(ρα)ρ′(rαβ) + 1 2ϕ ′_(rαβ ) )_rαβ i r αβ j rαβ ここで，各原子位置における原子応力を σα_ij = 1 V /N N ∑ β(̸=α) ( F′(ρα)ρ′(rαβ) + 1 2ϕ ′_(rαβ₎ )_rαβ i r αβ j rαβ (2.27) と定義すると，系の応力は σij = 1 N N ∑ α σ_ijα (2.28) となる．

2.5.4 原子弾性係数 弾性係数も応力と同様に U ≈ Etot の場合には Cijkl = 1 V ∂2Etot ∂ηij∂ηkl (2.29) であるので，平衡状態からの仮想均一変形を考えるとEAMポテンシャルにおける弾性係 数は以下のようになる． Cijkl = 1 2V N ∑ α N ∑ β(̸=α) ∂rαβ ∂ηkl ∂ ∂rαβ (∑N α N ∑ γ(̸=α) {F′_(ρα )ρ′(rαγ) + 1 2ϕ ′_(rαγ )}r αγ i r αγ j rαγ ) = 1 V [∑N α N ∑ β(̸=α) F′(ρα) ( ρ′′(rαβ)− ρ ′_(rαβ₎ rαβ ) r_iαβrαβ_j rαβ_k r_lαβ (rαβ₎2 + N ∑ α F′′(ρα) { ∑N β(̸=α) ρ′(rαβ)r αβ i r αβ j rαβ }{ ∑N γ(̸=α) ρ′(rαγ)r αγ k r αγ l rαγ } + 1 2 N ∑ α N ∑ β(̸=α) { ϕ′′(rαβ)− ϕ ′_(rαβ₎ rαβ }_rαβ i r αβ j r αβ k r αβ l (rαβ₎2 ] (2.30) 応力と同様に，各原子位置における原子弾性係数を以下のように定義する． C_ijklα = 1 V /N [ ∑N β(̸=α) F′(ρα) ( ρ′′(rαβ)− ρ ′_(rαβ₎ rαβ ) rαβ_i r_jαβrαβ_k r_lαβ (rαβ₎2 + F′′(ρα) { _∑N β(̸=α) ρ′(rαβ)r αβ i r αβ j rαβ }{ _∑N γ(̸=α) ρ′(rαγ)r αγ k r αγ l rαγ } + 1 2 N ∑ β(̸=α) { ϕ′′(rαβ)− ϕ ′_(rαβ₎ rαβ } r_iαβrαβ_j r_kαβrαβ_l (rαβ₎2 ] (2.31) これより，系の弾性係数は Cijkl = 1 N N ∑ α C_ijklα (2.32)

のように原子弾性係数の平均となる．

2.5.5 原子弾性剛性係数

以上で定義した原子応力，弾性係数から，原子弾性剛性係数は以下で評価できる．

B_ijklα = C_ijklα + (σα_ilδjk+ σjlαδik

+ σα_ikδjl+ σjkαδil− 2σijαδkl)/2 (2.33)

Wangらによる提案17) に従い，Voigt対称性をもたせたB_ijklα sym ≡ (B_ijklα + B_lkjiα )/2 を

用いて安定性評価を行う．なお，ここでのi∼lは直交座標の指標(1,2,3)もしくは(x,y,z)

であるが，以降ではxx，yy，zz，yz，zx，xyを1∼6とするVoigt表記を用いてBα_ij と表

3

モード

I

型き裂の引張シミュレーション

3.1

シミュレーション条件

図3.1に示すように，平板状の周期セル中央に貫通き裂を有する系を解析対象とする． き裂幅2cはセル寸法 Lx の1/5とした．図の3 つの結晶方位を対象とし，[001](010)， [¯110](001)，[¯1¯12](111)のき裂の解析を行った．き裂変位hはVoterのEAMポテンシャ ル15)でのNiのカットオフ半径0.479nmより大きく，かつき裂先端半径ρがグリフィス き裂の条件8a/π（aはNiの格子長さ0.352nm）を満たすよう，ρ=0.4nm，h=2ρ=0.8nm とする．き裂部の原子を取り除いた各モデルの原子総数はそれぞれ 178200, 179074, 173864である． 各モデルについて，全方向周期境界の下で垂直応力が0となるようセル辺長を調整しな がら100000fsの緩和計算を行った後，y 方向にひずみ増分∆εyy = 1.0× 10−6[/fs]を毎 ステップ与えることで引張シミュレーションを行う．数値積分の時間ステップは1fsであ るので，ひずみ速度に換算すると1.0× 109[/s]である．温度は熱搖動の影響を排除する ため0.1Kとし，速度スケーリングにより制御している．応力-ひずみ応答，ならびに，高 エネルギーの原子位置から決定したき裂長さ変化を求めると共に，1000fs毎(き裂進展近 傍では100fs毎)に記録した原子座標データから，各原子位置の原子応力 σα ij ならびに原 子弾性係数Cα ij からBijα を評価し，detBijα および固有値ηα(k)を求めた．なおdetBijα の 値は0[K]完全結晶における値8.8796× 1012[GPa]で無化元化している．

[100]

[001](010) crack

[010]

[001]

[-110]

[110]

[001]

[-110](001) crack

[111]

[-1-12](111) crack

[-1-12]

[-110]

3.2

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

3.2.1 応力-ひずみ線図と負のAES原子数変化 図3.2に応力ひずみ線図及び負のAES原子数の変化を示す．いずれのモデルにおいて も，き裂が存在するにも関わらず無負荷平衡状態ではdetB_ijα < 0の原子は存在しない． ひずみ_{εyy = 0.03}近傍からdetB_ijα が負となる原子が現れ始める．初めて負の原子が現れ た点を応力-ひずみ曲線上に×印で示した．また図中の○印は，可視化ソフトatomeye18)

のcentral symmetry parameterによりき裂先端からの欠陥生成（表面原子以外に生じた

欠陥）が確認された点を示している．応力-ひずみ線図を見ると，全てのモデルにおいて ○印の後から非線形性が現れており，非弾性変形の開始に対応していることがわかる． [001](010)model [-110](001)model [-1-12](111)model St re ss , σyy , [ G Pa ]

strain,

ε

_yy 0 0.1 0 5 10 15 [001](010)model [-110](001)model [-1-12](111)model

strain,

ε

_yy N um be r o f n eg at iv e de tB α ij a to m s 0 0.1 0 1000 2000

(i) stress-strain curves. (ii) Change in the number of detBα_ij < 0 atoms.

Fig.3.2 Stress-strain curves and change in the number of negative detBα

図3.3にシミュレーション中の各モデルの変形の様子をatomeye18)のMisesひずみで 色分けして示す．Misesひずみは0∼0.2の範囲で色分けしており，ひずみの値が大きい場 所ほど明るい色で表示されている．[001](010)モデルでは，応力-ひずみ線図で×で示した 点(図(a)(ii))をみると，き裂の先端に高いMisesひずみが確認される．その後，き裂先端 からの転位射出と共にき裂はぜい性的に進展した(図(a)(iii)，(iv))．ひずみεyy = 0.057 でき裂進展開始し，その後一定応力で破壊が進行する．[¯110](001)モデル，[¯1¯12](111)モ デルでも同様に非線形性が現れる点(図(ii))ではき裂先端に高いMisesひずみが観察さ れるが，その後は転位射出によってき裂が鈍化する(図(iii)).[¯110](001)モデルではき裂 から転位射出することで応力が鈍化し，εyy = 0.059でき裂進展開始する．[¯1¯12](111)モ デルではき裂の左側からは多数の転位が発生し，右側はき裂上端のみから転位が射出さ れた．

(i) ε_yy=0

(i) εyy=0

(i) εyy=0 (ii) εyy=0.030 (iii) εyy=0.051 (iv) εyy=0.072

(ii) εyy=0.041 (iii) εyy=0.059 (iv) εyy=0.086

(ii) εyy=0.048 (iii) εyy=0.057 (iv) εyy=0.072

(a)[001](010)crack (b)[110](001)crack (c)[112](111)crack

3.2.2 き裂先端の負のAES原子の分布 最初にdetB_ijα が負の原子が観察されたステップにおけるき裂近傍のスナップショット を図3.4に示す．図ではdetB_ijα < 0の原子を赤で，それ以外を緑で着色して示している． いずれのモデルでもき裂先端にdetB_ijα が負の原子が観察されており，奥行き方向(周期方 向)の原子列がdetB_ijα < 0となっている．[001](010)，[¯1¯12](111)モデルでは表面原子が 負となっているのに対して，[¯110](001)モデルでは表面から一原子内側の原子が負となっ

ている．また図3.5に，AtomeyeのCentral Symmetry Parameter(CSP)によって，き

裂表面以外は欠陥と判定される原子が存在しないステップでの負のAES原子の分布(左)

と，2000∼4000fs後で欠陥と判定された原子のスナップショット(右)を示す．ぜい性的

なき裂進展を示した[001](010)モデルでは，き裂先端に負のAES原子が多く密集してい

る．一方，転位射出によりき裂が鈍化した[¯110](001),[¯1¯12](111)き裂では負のAES原子

を起点として転位が射出されている．

(i)[001](010) model (ii)[110](001) model- (iii)[112](111) model-

-Fig.3.4 Snapshots of atoms at emerging point of negative AES atoms.

図3.6 は応力-ひずみに非線形性が現れてしばらく後のスナップショットを，AES と

その前縁に負のAES原子が多くみられる．[¯110](001)き裂では転位射出後に積層欠陥は 消滅するので，先の図3.5のすべり面に沿った欠陥はなくなり，き裂底に生じたステップ 部分にのみCSPで欠陥が認められる．これに対応して，負のAES原子は図3.5のように き裂の底部だけではなくステップ中央にも生じている．[¯1¯12](111)き裂では，き裂右側か ら射出された完全転位(右上)が不安定となっているが，左側ではき裂先端から発生した 部分転位のループ(積層欠陥部)が途中で止まっている．この転位の前縁には負のAES原 子がみられないため成長が止まった可能性がある．

(a) [001](010) crack (b) [110](010) crack (c) [112](111) crack -(i)ε_yy=0.048 (ii)ε_yy=0.052 (i)ε_yy=0.043 (ii)ε_yy=0.045 (i)ε_yy=0.029 (ii)ε_yy=0.031

(a) ε=0.058([001](010) crack) (b) ε=0.06([110](001) crack) (c) ε=0.42([112](111)crack)

-Fig.3.6 Snapshots of negative AES atoms(left) and defect atoms visualized by central cymmetry parameter(right).

3.2.3 固有値の変化 B_ijα_{εj = η}α_εiの固有方程式の解 ηα(k) (k = 1∼6)を各原子について求め，その変化 を調べた．図3.7に引張シミュレーション中の第一固有値が負の原子数の変化を示す．初 めて第一固有値が負の原子が現れる点はdetB_ijα < 0の原子が現れる点と完全に一致して いる．また図3.2(ii)に示したdetBα_ij < 0の原子数変化のグラフとほぼ同じ概形となって いることから，不安定原子の大部分は第一固有値のみが負の原子を示していることがわか る．図3.8は，各ひずみでηα(1)_{が最も小さい原子を調べ，その値をプロットしたものであ} る．図中にはき裂長さの変化を青線で示してある．また図中の×印は初めてdetB_ijα < 0 の原子が現れる点である．いずれも×印の前から減少しており，[001](010)，[¯1¯12](111) き裂ではほぼ-40GPa(単位はB_ijα_{εj = η}α_{εi = σ}i より)で一定の値をとった後，き裂の長 さが変化している．[¯110](001)き裂では一度正値に回復するような複雑な挙動を示してい [001](010)model [-110](001)model [-1-12](111)model

strain,

ε

_yy

N

um

be

r o

f n

eg

at

iv

e

η

α (1 )

at

om

s

0 0.1 0 1000 2000

るが，転位射出によってき裂近傍の不安定性が回復した可能性がある．き裂長さが増加す る際はηα(1) の最小値は約半分の-20GPa程度となっている．表3.1に×印，非弾性変形 開始点(CSPにより欠陥原子を確認した点)におけるひずみの値と，そのときの負のηα(1) を有する原子数，最小のηα(1) の値をまとめて示した．×印における負の原子数，ηα(1) の値に大きな差は無く，局所不安定が発生するときの力学条件はモデルによる違いがない と考えられる．一方，その後欠陥となって実際に構造変化を生じる点では負の原子数(領 域)には大きな違いがあるが，ηα(1)の最小値(変形モードに対応)は-15∼-35GPaの範囲 内にある． strain, εyy M in im um o f e ig en va lu e, η α (1 ) C rac k l en gth Minimum of eigenvalue Crack length 0 0.05 0.1 -60 -40 -20 0 20 0.2 0.4 0.6 strain, εyy M in im um o f e ig en va lu e, η α (1 ) C rac k l en gth Minimum of eigenvalue Crack length 0 0.05 0.1 -60 -40 -20 0 20 0.2 0.4 0.6 strain, εyy M in im um o f e ig en va lu e, η α (1 ) C rac k l en gth Minimum of eigenvalue Crack length 0 0.05 0.1 -60 -40 -20 0 20 0.2 0.4 0.6

(a)[001](010) crack (b)[110](001) crack

(c)[112](111) crack

Table 3.1 Strain, number of ηα(1) < 0 atoms, and minimum of ηα(1) at emerging point of negative AES and defect.

emergence point of ηα(1) < 0 defect nucleation by CSP

εyy atoms min. of ηα(1) εyy atoms min. of ηα(1)

[001](010) crack model 0.035 40 -5.08 0.050 294 -35.91 [¯110](001) crack model 0.042 56 -2.99 0.043 168 -14.63 [¯1¯12](111) crack model 0.026 48 -3.31 0.030 69 -22.44 図3.9は引張シミュレーション中に第二固有値が負となった原子数変化を示したもので ある．図中にはき裂長さの変化を青線で示してある．[001](010)き裂では，第二固有値が 負の原子数がパルス状に増加した直後にき裂が進展した．[¯110](001)モデルでは非弾性変 形の開始点(εyy = 0.043) では第二固有値が負の原子は観察されず，εyy = 0.07近傍で 初めて第二固有値が負の原子が現れる．青線のき裂長さはこの点から急激に増加してい る．[¯1¯12](111)き裂では他に比べ早くから第二固有値が負の原子がみられるがその数は少 ない．転位射出による鈍化のため，き裂長さ変化は他に比べ緩やかであるが第二固有値が 負の原子が出現した点からき裂の成長が始まっており，き裂の進展は第二固有値が負の原 子によって測定できることが示唆される．

0 0.05 0.1 0 10 20 30 40 50 0.2 0.4 0.6 0 0.05 0.1 0 10 20 30 40 50 0.2 0.4 0.6 Number of negative η α(2) atoms Number of negative η α(2) atoms Number of negative η α(2) atoms

crack length crack length

crack length (i)[001](010)model (ii)[110](001)model (iii)[112](111)model -0 0.05 0.1 0 10 20 30 40 50 0.2 0.4 0.6

図 3.10に第二固有値が負の原子が出現したステップでのき裂近傍のスナップショッ トを，ηα(1) < 0の原子を赤，ηα(2) < 0の原子を青，それ以外を緑で示す．第二固有値 が負となった原子は第一固有値も負となっている．したがって，detBα_ij の判定では正と 判定される原子であることに注意が必要である．[001](010)き裂ではき裂先端の表面側 に，第二固有値が負の原子が二列(奥行き方向全ての原子)，き裂面上下に発生している． [¯110](001)き裂では転位射出により生じた図中右下の入り込み部に一列，また第一固有 値が負の原子列が凝集した右上結晶内部に一列，第二固有値が負の原子列が認められる． [¯1¯12](111)モデルでは同様に転位射出による右上の入り込み部に第二固有値が負の原子列 があり，この部分はわずかにき裂進展した．左下にも第二固有値が負の原子がみられる が，この部分では奥行き方向に連続しておらずき裂進展は生じていない．

(i)ε=0.048 (ii)ε=0.071 (iii)ε=0.061

4

ダブルき裂モデルでのシミュレーション

4.1

シミュレーション条件

図4.1に示すように，前章で用いたセルの四隅に4分の1貫通き裂を配置したモデル で解析を行った．き裂幅2cは前章と同様にセル寸法 Lx の1/5とした．き裂部の原子を 取り除いた各モデルの原子総数はそれぞれ178200, 177828, 173008である．前章と同様 に，全方向周期境界の元で垂直応力が0となるようセル辺長を調整しながら温度0.1Kで 100000fsの緩和計算を行った後，y 方向にひずみ増分∆εyy = 1.0× 10−6[/fs]を毎ステッ プ与えることで引張シミュレーションを行う．その後の引張条件についても前章と同様と した．

[100]

[001](010)crack

[010]

[001]

[-110]

[110]

[001]

[-110](001)crack

[111]

[-1-12](111)crack

[-1-12]

[-110]

4.2

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

4.2.1 応力-ひずみ線図と負のAES原子数変化 図4.2に応力ひずみ線図及び負のAES原子数の変化を示す．初めて負の原子が現れた 点，ならびにCSPで欠陥が認められた点を応力-ひずみ上にそれぞれ×，○印で示した． [110](001)き裂，[112](111)き裂は前章と同様，○印以降非線形性が現れた後応力上昇が 鈍化し，弾完全塑性体のような挙動を示している．一方，[001](010)き裂では応力-ひずみ 曲線はピーク後著しく低下した．図4.3に変形の様子を，負のAES原子を赤，それ以外 を緑で着色して示す．図(i)は各モデルにおける応力-ひずみ線図の×印近傍でのスナップ ショットを示している．その後[001](010)き裂では，き裂先端にAESが負の原子が集中 [001](010)model [-110](001)model [-1-12](111)model

strain,

ε

_yy N um be r o f n eg at iv e de tB α aij to m s 0 0.1 0 1000 2000 [001](010)model [-110](001)model [-1-12](111)model St re ss , σyy , [ G Pa ]

strain,

ε

_yy 0 0.1 0 5 10

(i) stress-strain curves. (ii) Change in the number of detBα_ij < 0 atoms.

Fig.4.2 Stress-strain curves and change in the number of negative detBα

して生じた後に，四隅のき裂がぜい性的に成長している(図(iii),(iv))．変形後期では四隅 のき裂が合体して材料がほぼ分離した状態となったため先の応力急減を生じたものと結論 づけられる．[¯110](001)，[¯1¯12](111)き裂ではき裂先端に負のAES原子が現れるが，前章 と同様に転位射出によるき裂鈍化を示した．(図(iv))． (i)ε=0.037 (ii)ε=0.051 (i)ε=0.045 (ii)ε=0.067 (a)[001](010)crack (b)[110](001)crack (i)ε=0.027 (iii)ε=0.054 (c)[112](111)crack -(iii)ε=0.070 (iv)ε=0.097 (iv)ε=0.086 (iv)ε=0.113 (ii)ε=0.040 (iii)ε=0.077

4.2.2 負のAES原子が発生するひずみと分布形態 負の AES原子が発生したステップでのひずみの値を表4.4にまとめた．いずれのモデ ルもわずかに単一き裂の場合のほうがひずみが小さいが，ほぼ同じ値となっている．複数 き裂の方が応力遮蔽によって局所不安定に達するのが単一き裂より遅くなったものと考え られる．負のAES原子発生時の中央のき裂と四隅のき裂周囲での負の原子分布を，前章 で用いた単一き裂モデルと合わせて図4.4に示す．図中では負の原子を赤，それ以外を緑 で示している．[001](010)き裂では，単一き裂時にはき裂の両端に奥行き方向(周期方向) に整列して負のAES原子が発生したが，複数き裂の場合では中央のき裂の右下と，四隅 のき裂の左下にのみ観察され，また奥行き方向にも全ての原子列が負になってはいなかっ た．ただし，この直後のステップ(1000fs後)のスナップショットでは複数き裂でも単一 き裂と同様の分布となった．[¯110](001)き裂では，単一き裂では一原子内部の原子のみが 負となったのに対して，複数き裂ではどちらのき裂も表面原子と，その一つ内部の原子が 負のAESを示し，単一き裂時に負の値を示した原子列と異なる原子列が負となっている． [¯1¯12](111)き裂では，単一き裂ではき裂の先端原子とその上下の原子列が負となったが複 数き裂では先端原子のみが負となった．ただし，[001](010)き裂の場合と同様，1000fs後 の原子配置では負の原子分布は単一き裂と同じになった．

Table 4.1 Strain at emerging point of negative AES atoms.

single crack model double crack model

[001](010) 0.035 0.037

[¯110](001) 0.042 0.044

(a)single crack (b)double crack (left:center right:corner) (i)[001](010)crack (ii)[110](001)crack (iii)[112](111)crack

表4.2に各モデルでの非弾性変形開始直前 (CSPでの欠陥生成点，非線形性が現れる 点)のひずみと負のAES原子数を，前章で用いた単一き裂のものと合わせて示す．ただし [¯110](001)き裂ではき裂の右端と左端で欠陥原子が出るタイミングが異なるため，それぞ れ欠陥原子が出る直前のステップでの原子数を求め，足したものを用いている．鉄による 先行研究9)では，ほとんど塑性変形を生じずへき開する結晶方位([001](110),[¯1¯12](111)) では，き裂の大きさや数，配置に関わらずき裂先端に現れる負のAES原子数が一定であ ることが明らかにされている．fcc-Niでは[001](010)き裂では非弾性変形時にき裂が進 展し，[¯110](001)，[¯1¯12](111)き裂では転位射出によるき裂の鈍化が起こった．しかしい ずれのモデルにおいても非弾性変形が起こるひずみがほぼ同じで，ダブルき裂モデルは単 一き裂モデルのほぼ二倍の負のAES原子を生じている．鉄の場合と同様に，き裂先端の 不安定領域は一定の寸法である可能性がある．

Table 4.2 Number of negative AES atoms at the onset of inelastic deformation.

single crack model double crack model

strain atoms strain atoms

[001](010) 0.050 294 0.053 602

[¯110](001) 0.044 168 0.045 266

4.2.3 固有値による検討 図4.5に最小固有値(ηα(1)のうち最も小さいもの)の変化を示す．図中にはき裂の長さ の変化を赤線と青線で示してある．[001](010)き裂では_{εyy = 0.025}近傍から最小固有値 の値が前駆的に減少し始める．その後すぐにき裂長さにステップ上の変化があらわれてい るが，このときはまだηα(1)は負ではない．ひずみεyy = 0.025以降ηα(1)は負値となる が，ηα(1)が-40GPa(B_ijα_{εj = σ}α_i = ηα_εiより)程度になったとき，四隅のき裂が急激に成 長している．その値を保った状態となるが，この点からき裂の長さが急激に伸び始めてい る．．[¯110](001)き裂ではやはり×印の少し前から最小値が減少し始める．この結晶方位 ではき裂開口せず転位射出するが，負値に達した後のピーク状の変化が転位の射出に対応 するものと考えられる．[¯1¯12](111)き裂ではηα(1) の減少が最も急激で，き裂成長時-40 ∼-60GPaの一定の値を示す．図4.6に第二固有値が負の原子数の変化を，き裂長さの変 化と共に示す．前章と同様，[001](010)き裂，[¯110](001)き裂では負の原子数がパルス状 の増加を示したのちにき裂長さが変化している．[¯1¯12](111)き裂も他より早い段階で負の 原子が観察されるとともに，穏やかにき裂長さが変化している．ダブルき裂モデルにおい てもき裂長さが変化し始める点と第二固有値が負の原子が発生する点が対応している．

(a)[001](010)crack (b)[110](001)crack (c)[112](111)crack -strain,

ε

_yy M in im um o f e ig en va lu e, η α (1 ) minimum of eigenvalue cracklength(center crack) cracklength(corner crack) C rac k l en gth 0 0.05 0.1 -60 -40 -20 0 20 0.2 0.4 0.6 strain,

ε

_yy M in im um o f e ig en va lu e, η α (1 ) minimum of eigenvalue cracklength(center crack) cracklength(corner crack) C rac k l en gth 0 0.05 0.1 -60 -40 -20 0 20 0.2 0.4 0.6 strain,

ε

_yy M in im um o f e ig en va lu e, η α (1 ) minimum of eigenvalue cracklength(center crack) cracklength(corner crack) C rac k l en gth 0 0.05 0.1 -60 -40 -20 0 20 0.2 0.4 0.6

strain, ε_yy N um be r o f n eg at iv e η α (2 ) at om s C rac k l en gth minimum of eigenvalue cracklength(center crack) cracklength(corner crack) 0 0.05 0.1 0 25 50 75 0.2 0.4 0.6 strain, ε_yy N um be r o f n eg at iv e η α (2 ) at om s C rac k l en gth minimum of eigenvalue cracklength(center crack) cracklength(corner crack) 0 0.05 0.1 0 25 50 75 0.2 0.4 0.6 strain, ε_yy N um be r o f n eg at iv e η α (2 ) at om s C rac k l en gth minimum of eigenvalue cracklength(center crack) cracklength(corner crack) 0 0.05 0.1 0 25 50 75 0.2 0.4 0.6

(a) [001](010) crack (b) [110](001) crack

(c) [112](111) crack

5

厚さ方向自由度の影響

5.1

シミュレーション条件

本章では第3章で用いたモデルを厚さ方向に1.5倍して，き裂前縁方向の自由度を増や した系を解析対象とする．平板状の周期セル中央に貫通き裂を導入し，き裂寸法，結晶方 位は第3 章のものと同様の条件である．き裂部の原子を取り除いた原子総数はそれぞれ 249480, 268611, 260796である．各モデルについて，全方向周期境界の下で垂直応力が0 となるようセル辺長を調整しながら100000fsの緩和計算を行った後，y 方向にひずみ増 分∆εyy = 1.0× 10−6[/fs]を毎ステップ与えることで引張シミュレーションを行う．

5.2

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

5.2.1 応力-ひずみ線図と変形の様子 図5.1に応力ひずみ線図及び負のAES原子数の変化を示す．図中下には比較のため第 3章の図を再掲する．いずれのモデルにおいても無負荷平衡状態では負の原子は存在して おらず，ε = 0.03近傍からdetB_ijα が負の原子が現れ始める．初めて負の原子が現れた点 を応力-ひずみ曲線上に×印で，CSPで欠陥を確認したひずみを○で示す．○印から応力 -ひずみの非線形性が現れるという点は第3章，第4章の結果と同様である．[¯1¯12](111)き 裂の応答が，第3章のそれよりわずかに変化しているが，それ以外はほぼ同じである．

[001](010)model [-110](001)model [-1-12](111)model St re ss , σyy , [ G Pa ]

strain,

ε

_yy 0 0.1 0 5 10 15 [001](010)model [-110](001)model [-1-12](111)model

strain,

ε

_yy N um be r o f n eg at iv e de tB α ij a to m s 0 0.1 0 1000 2000 3000

(a) stress-strain curves. (b) Change in the number of detB_ijα < 0 atoms.

thick model [001](010)model [-110](001)model [-1-12](111)model St re ss , σyy , [ G Pa ]

strain,

ε

_yy 0 0.1 0 5 10 15 [001](010)model [-110](001)model [-1-12](111)model

strain,

ε

_yy N um be r o f n eg at iv e de tB α ij a to m s 0 0.1 0 1000 2000

(c) stress-strain curves. (d) Change in the number of detB_ijα < 0 atoms.

base model

図5.2に引張シミュレーション中の変形の様子を，負のAES原子を赤，それ以外を緑で 着色して示す．[001](010)き裂のみぜい性的にき裂が進展し，[¯110](001)，[¯1¯12](111)き 裂ではき裂先端に負の原子が現れた後転位が射出され転位が鈍化するという点は第3章と 同じである．図5.2(i)はそれぞれのモデルで最初に負のAES原子が現れたステップでの スナップショットを示しているが，その分布は図3.4と一致している．しかし[¯1¯12](111) き裂では，第3章ではき裂左端の負のAES原子は転位ループが途中で止まっていたが， 厚さ方向が大きい今回の結果ではき裂左端からも完全転位として負のAES原子が移動す る様子が観察された．図5.3に[¯1¯12](111)き裂での，図5.1の○印の前のステップ(表面 原子以外欠陥がない)における負のAES原子の分布(左)と，その後欠陥と判定された原 子のスナップショット(右)を示す．[¯1¯12](111)き裂ではき裂前縁長さ方向の自由度を増 やしたことによって，両方のき裂先端から完全転位が射出されている．

(i) ε_yy=0.035 (ii) ε_yy=0.048 (iii) ε_yy=0.062

(a)[001](010)crack

(b)[110](001)crack

(i) ε_yy=0.043 (ii) ε_yy=0.060 (iii) ε_yy=0.070

-(i) ε_yy=0.026 (ii) ε_yy=0.035 (iii) ε_yy=0.049

(c)[112](111)crack-

(a)ε=0.03 (b)ε=0.034

Fig.5.3 Snapshots of [¯1¯12](111) crack negative AES atoms(left) and defect atoms

visualized by central cymmetry parameter(right).

5.2.2 負のAES原子数と固有値 図 5.1 の×印におけるひずみと負の原子数を，第 3 章における値と共に表 5.1 に示 す．いずれのモデルも負の原子が初めて現れるひずみはほぼ等しく，またその数も1.5倍 ([001](010)き裂のみ1.4倍)となっていることから，detBα_ij < 0の原子が発生するときの 力学条件には厚さ方向自由度の影響はない．図5.4に最小のηα(1)の変化を，第3章のも のと併せて示す．両者はほとんど同じで，[001](010)き裂，[¯1¯12](111)き裂では-40GPa で一定の値をとり，[¯110](001)き裂では-20GPa程度の値を示して変形が進行する．また ηα(2) _{が負になった原子数の変化を，第3章のものと図5.4併せて示す．[001](010)き裂} ではパルス状の増加を示した際の原子数が増えているが，それ以外はほぼ同じである． [¯110](001)き裂はηα(2) < 0の原子が出現するひずみがわずかに早くなっており，また2 番目のパルスが非常に大きくなっている．[¯1¯12](111)き裂では，第3章ではεyy = 0.03 近傍からηα(2) < 0の原子が現れていたが，それが現れておらず変形後期にパルス状に現

[001](010)model [-110](001)model [-1-12](111)model M in im um o f e ig en va lu e, η α (1 )

strain,

ε

_yy 0 0.02 0.04 0.06 0.08 -40 -20 0 20 [001](010)model [-110](001)model [-1-12](111)model

strain,

ε

_yy M in im um o f e ig en va lu e, η α (1 ) 0 0.02 0.04 0.06 0.08 -40 -20 0 20

(a) base model

(b) thick model

Fig. 5.4 Change in the minimum 1st eigenvalue.

度によって変形モードが異なる可能性がある．各モデルでのηα(2) < 0 の原子数の変化

を，き裂長さ変化と併せて図5.6に示す．[001](010)，[¯110](001)き裂は ηα(2) _{< 0の出}

現とき裂長さが急増する点が対応しているが，[¯1¯12](111)き裂はき裂長さ変化が緩慢であ

るためηα(2)< 0の出現と対応していると結論づけることは難しい．

Table 5.1 Strain at emerging point of negative AES atoms.

single crack model thick model

strain atoms strain atoms

[001](010) 0.035 40 0.035 56

[¯110](001) 0.042 56 0.043 84

(a) base model

(b) thick model

[001](010)model [-110](001)model [-1-12](111)model

strain,

ε

_yy

N

um

be

r o

f n

eg

at

iv

e

η

α (2 )

at

om

s

0 0.02 0.04 0.06 0.08 0 20 40 60 80 [001](010)model [-110](001)model [-1-12](111)model

strain,

ε

_yy

N

um

be

r o

f n

eg

at

iv

e

η

α (2 )

at

om

s

0 0.02 0.04 0.06 0.08 0 20 40 60 80

Number of negative ηα(2) _atoms Crack length N um be r o f n eg at iv e η α (2 ) at om s C rac k l en gth strain, εyy 0 0.05 0.1 0 20 40 60 80 0.2 0.4 0.6

Number of negative ηα(2) _atoms

Crack length N um be r o f n eg at iv e η α (2 ) at om s C rac k l en gth strain, εyy 0 0.05 0.1 0 20 40 60 80 0.2 0.4 0.6

Number of negative ηα(2) _atoms

Crack length N um be r o f n eg at iv e η α (2 ) at om s C rac k l en gth strain, εyy 0 0.05 0.1 0 20 40 60 80 0.2 0.4 0.6 (a)[001](010)crack (b)[110](001)crack (c)[112](111)crack

6

結論

本研究では，Niの薄板周期セルにモードI型き裂を導入したモデルに，分子動力学法に より引張シミュレーションを行い，き裂進展挙動を原子弾性剛性係数B_ijα の行列式detB_ijα と固有値ηα の観点から議論した．以下に本論文で得られた結果を総括する． 第2章では解析手法の概要を述べた．はじめに分子動力学法の基礎方程式を示し，原子 間相互作用の評価に用いたEAMポテンシャルの概要を説明した．さらに，各原子の安定 性評価の指標となる原子弾性剛性係数について概説した． 第3章では薄板状の周期セル中央に貫通き裂を有する系を解析対象とし，[001](010)， [¯110](001)，[¯1¯12](111)の三種類の結晶方位について検討した．いずれのモデルにおいて も無負荷平衡状態ではdetB_ijα < 0の原子は存在しておらず，き裂近傍に detB_ijα < 0と なる原子が生じた後，き裂先端から欠陥が発生し応力-ひずみ線図に非線形性が現れた． detB_ijα < 0の原子が初めて出現するときのひずみはそれぞれ0.035，0.042，0.026で，き 裂先端近傍の原子列がそれぞれ2，2，3本，detBα ij < 0となっている．固有値の変化を 調べた結果，構造変化が起こる点(欠陥が発生する点)ではηα(1)の最小値がほぼ一定の負 の値([001](010)，[¯1¯12](111)は-40GPa，[¯110](001)は-20GPa）となっていた． 第4章では，先の3つの周期セルの四隅に4分の 1き裂を配置したダブルき裂モデル のシミュレーションを行い，周期き裂の配置を変えた影響を調べた．初めてdetB_ijα < 0 の原子が出現するひずみは先の単一き裂に比べてわずかに高ひずみ側(0.002程度)となっ たが，これはき裂配置による応力遮蔽効果と考えられる．欠陥生成時のdetBα ij < 0原子 の数はいずれも単一き裂モデルのほぼ二倍であり，き裂先端の不安定領域は一定の寸法で

ほぼ先述の値(-40GPaと-20GPa)と一致していた． 第5章では第3章のシミュレーションセルを厚さ方向に1.5倍とし，き裂前縁方向の自 由度の影響について検討した．[¯1¯12](111)き裂のみ，厚さ方向の自由度を増やしたことに よって変形モードがわずかに変化することが，応力-ひずみ，原子配置変化，第二固有値の 変化などから示された．ただしdetBα ij < 0の原子が初めて現れるひずみにはどのモデル も変化がなく，その原子数も1.5倍([001](010)のみ1.4倍)となっていることから，き裂 先端にdetB_ijα < 0の原子が発生するときの力学条件は厚さ方向の自由度の影響はない． 欠陥生成してき裂が開口するときの第一固有値の値も先述の値と同じである．以上を総合 すると，結晶方位による差はあるが，き裂先端に初めてdetB_ijα < 0 が出現するときの局 所力学状態は常に同じ(detB_ijα < 0の原子列＝不安定領域が一定)であり，またその後き 裂先端から転位射出しき裂が鈍化・開口するときのBα ij の第一固有値が一定の値を示すこ とが本研究で得られた事実である．

参考文献

[1] 屋代・冨田，日本機械学会論文集(A編),67，656,(2001)，678

[2] 屋代・西村・冨田，材料，54,10,(2005)，1053

[3] 屋代・藤原，材料，60,11,(2011),968

[4] Yashiro, K. and Yamane, T., Proceedings of ICCM2014, ISSN2374- 3948(online), Sientech Publisher (2014).

[5] Chang R, Int. Journ. of Fracture Mech., 6,(1970),pp111-125 [6] Gumbsch P, Mater. Sci. Eng.:A,Vol.319-321,(2001),pp.1-7

[7] Shimada, T., Ouchi, K., Chihara, Y. and Kitamura, T., Scientific Reports, No.5(2015),p.8596.

[8] Yashiro, K., Tsugawa Y., and Katayama H., Advanced Structured Materials, From Creep Damage Mechanics to Homogenization Methods, Springer, pp.557-570(2015)

[9] 屋代・片山・西川・吉原，日本機械学会論文集, Vol.81, No.829, 15-00271(2015)

[10] 上田顯, コンピューターシミュレーション, (1990), 朝倉書店.

[11] 上田顯, 分子シミュレーション-古典系から量子系手法まで-, (2003), 裳華房.

[12] Daw, M.S., and Baskes, M.I., Physical Review Letters, 50,(1983),1285 [13] Daw,M.S., and Baskes,M.I., Physical Review B, 29,(1984),6443

[14] Parisi, G., and Wu, Y., The Science of Sin, 24, 483, (1981)

[15] Voter, A.F., Chen, S.P., Siegel, R.W., Weertman, J.R., Sinclair, R. (Eds.), Char-acterization of Defects in Materials, 82,(1987),175

[17] Wang, J., Tip, S., Phillpot, S.R., and Wolf, D.,Phys.Rev.lett., 71-25, (1993), 4182.

謝辞

本研究を遂行するにあたり，仔細に渡る懇切丁寧な御指導を賜りました屋代如月教授 (現 岐阜大学)に，心より深く感謝の意を表します．また，本研究に対して，貴重な御助 言を賜りました田中克志教授，長谷部忠司准教授に，日々の研究活動において数々の励ま しの言葉を頂いた片山雷太技術職員に深く感謝致します． 本研究室配属後，若輩な私に沢山のアドバイスを下さり，基礎的な事まで懇切丁寧に教 えてくださった表面・界面工学研究室の諸先輩方に厚く御礼申し上げます． また，互いに切磋琢磨しながら研究に取り組み，苦楽を共にした同輩の穴瀬正典君，池 田翔伍君，岩中拓夢君，大橋正博君，佐藤博紀君，向井良多君に御礼申し上げます． 最後に，6年間の大学生活をあらゆる面から支えて頂いた家族には，深く感謝致します． ありがとうございました． 平成28年2月 片山 寛