強変形による結晶粒微細化過程の分子動力学解析と高効率プロセスの開発

研究代表者

九州工業大学情報工学部助手

松本龍介

１． はじめに多結晶体の結晶粒を微細化することで強度が上昇することは，以前からHall‑Petchの関係

1, 2)として知られてきた．近年になって，結晶粒サイズが1µm以下の微細粒材料を創成できる方法として，アモルファス金属からの結晶化を用いる方法3)，メカニカルアロイング（MA）法4)，強変形（Severe Plastic Deformation; SPD）を用いる方法5)等が開発されており，多くの材料において降伏応力がGPaのオーダーになることが示されている．本研究では，合金化を必要としないために材料のリサイクル性に優れ，さらに，比較的大きな寸法の材料の作成が可能である強変形を用いる方法に着目する．これまでに，強変形を用いて微細粒材料を創成する手法として，ECAP（Equal Channel Angular Pressing）6)，HPT（High PressureTorsion）7)，ARB（Accumulative Roll Bonding）8)等の様々なプロセスが提案されている．
しかしながら，実験的に強変形加工中に生じる内部組織変化を直接的に観察することは極めて困難であり，強変形によって結晶粒が微細化されるメカニズムの詳細は十分に明らかにされていない．そのため変形方向や温度等の加工条件の設定は，試行錯誤に頼っているところが大きい．そこで，計算機シミュレーションによる現象の解明に期待が寄せられている．分子動力学法（Molecular Dynamics; MD）9,10,11)は，計算機を用いて個々の原子に対する運動方程式を解くことで，系の時間発展を原子レベルで追跡する手法であり，それによって観察がほぼ不可能な強変形による結晶粒微細化過程の詳細を直接的に解析することが可能である．本研究では，MDシミュレーションを用いて，最も代表的な強変形プロセスの一つであるECAPによって生じる変形と等価な強せん断変形を取り扱い，結晶粒微細化の基本メカニズムを明らかにし，効率的な強変形プロセスを開発することを目的とする．
２． 解析方法と解析条件
2.1. ECAP の概要
ECAPでは図1(a)に示すような屈曲した金型に棒材を押し込み材料全体に単純強せん断変形を加えることで結晶粒を微細化する．図中で，座標系O‑XYZは金型に固定されている．一方，座標系o‑xyzは被加工材を基準として定義されており，せん断変形によって歪むことはなく，押し込み加工中は，材料挿入時の方向に固定されているものとする．押し込み加工によって，金型の出口から出てきた材料を図1(b)のように元の向きに戻し，つまりY軸まわりにΦ回転し，さらに，X，Y，Z軸まわりにそれぞれθX，θY，θZ回転して再び金型に押し込むことで加える単純強せん断の方向を変えることができる．押し込み加工ごとに，試料に固定された座標系o‑xyzは，座標系O‑XYZに対し，X軸まわりにθX，Y軸まわりに(Φ+θY)，Z軸まわりにθZ回転する．ここでは，最も一般的なΦ=90°のL字型の金型を用いたECAPを取り扱う．
この場合は一度の押し込み加工で工学ひずみγ=2が加わる．なお，実際に被加工材料をX軸，及び，Y軸まわりに回転するためには，形状に制限が加わり困難であるが，MDシミュレーシ
ョンでは実現可能であり，内部組織の応答を一般的に用いられている加工ルートと比較することで，より高効率なプロセスを模索することが可能である．

Fig.1. Schematic diagram of ECAP
2.2. 原子間ポテンシャル
本研究では，被加工材料として，α鉄を取り扱うこととする．また，bcc材料は，一般的にfcc材料に比べて内部組織変化が複雑であり，原子構造論的アプローチが有効であると考えられるが，これまでに結晶粒微細化過程のMD解析はほとんど行なわれていないことからも有意義である．本研究では，多原子かつ長時間の計算が必要であることから，精度と計算時間のバランスに優れるEAMタイプのポテンシャルを用いる．ここでは，SimonelliらがフィッティングしたパラメータセットA12)を用いる．本EAMポテンシャルでは，パラメータは凝集エネルギー，空孔形成エネルギー，弾性定数，及び，格子定数が実験値に一致するように設定されている．
2.3.初期構造の作成
本研究では，図1(a)中で灰色で示した材料内部の正方形領域をモデル化して取り扱う．解析対象として，板厚（Ly）を十分に小さく取り，計算セルの全方向に周期条件を適用した準3次元の多結晶体モデルを用いる．これは，金型がY軸方向に十分に長く，さらに試料内でY軸方向の変形がほとんど生じないという条件に近い．以下にその作成手順を説明する．
１． 結晶中心位置の決定（図2(a)を参照）

x軸方向寸法とz軸方向寸法がL=76.45nmの正方形領域を考え，その領域をx方向にnx個，z方向にnz個，つまりnx×nz個の長方形領域に分割し，各領域の頂点及び各領域の中心に初期結晶中心を配置する．次に各初期結晶中心を半径e×hの領域（鎖線で描かれた円内）で乱数を用いてランダムに動かし，結晶中心位置を決定する．ここで，h は初期結晶中心の最小距離である．本研究ではnx=4，nz=2とし，結晶粒数はn=nx×nz×2=16個，h=L/nxである．nxとnzは，図中で点線で示したように，1個の初期結晶中心を取り囲む6個の初期結晶中心の配置が正六角形に近くなるように選んでいる．また，結晶中心の分散状態を決定するパラメータeとして0.5を用いた．
前ステップで決定した結晶中心位置を代表点としてVoronoi分割を２． 結晶粒の作成することで個々の結晶粒が占める空間を決定し，それぞれの領域にy軸方向に<1,1,0>をとり，乱数を用いてy軸まわりに回転した完全結晶をはめこむことで多結晶体を作成する．y軸方向に<1,1,0>を取ったことで，xz面内にすべり方向である<1,1,1>が存在することになる．ここでは，板厚（y軸）方向には，10層の原子を並べている．
３． 近接原子の消去

人為的に結晶粒を並べたため，粒界近傍で原子間距離が非常に近い原子の組み合わせが存在する．そこで，最隣接原子間距離のf=0.73倍以下の原子の組み合わせに対して，片方の原子をランダムに消去する操作を行う．
４． 緩和計算

以上の手順により作成された多結晶体モデルに対して，垂直応力がσx=σy=σz=‑0.1GPaになるように計算セルのサイズを制御しつつ，300Kの一定温度条件下で5psの緩和計算を実施する．緩和計算後の辺寸法，平均粒径，原子数を表1に示す．なお，本研究で用いるMD計算では，時間積分法には，速度Verlet法を用い時間ステップはΔt =2fsとする．また，温度制御法には速度スケーリング法を用い，高速化・大規模化のためにブロック分割法とリンクセル法を用いる．最終的に得られた初期多結晶構造を図2(b)に示す．図中では，bcc構造以外の構造を有する原子（粒界部に存在）を黒色で示し，bcc構造を有する原子に対しては結晶粒を区別しやすいように結晶方位に応じて色を付けて示している．

Fig.2. Analysis model2.4. Lees‑Edwards周期境界条件
Φ=90°のECAPによる結晶粒微細化過程の解析を実施するためには，上記のように工学ひずみγ=2の単純せん断変形を繰り返し加える必要がある．しかしながら，周期境界条件下でよく用いられる計算セルの形状を変形させる方法では，大変形時にブロック分割法の効率が下がる等の不都合を生じる．そこで，ここではLees‑Edwards（L‑E）周期境界条件により，ECAPで加わる単純強せん断変形と等価な変形を実現する．L‑E周期境界条件を用いて，例えば，x軸方向に単純せん断変形を加える場合には，計算セルの上下（±z軸方向隣）のレプリカセルをそれぞれx軸方向と‑x軸方向に速度uで動かす．つまり，せん断変形は，z軸方向の周期境界をまたぐ相互作用や原子運動に対して，セルの変位Δx を考慮することで実現できる．さらに，毎ステップ，計算セル内の原子にひずみ速度に対応する変位勾配を与えることで，境界での変位の不連続をなくすことができる．また，L‑E周期境界条件では，上下のレプリカセルの相対的な移動量が，計算セルのせん断方向の寸法（Lx）の整数倍である
場合は，x軸方向の単純せん断変形の後に，左右（±x軸方向隣）のレプリカセルを上下に動かすことが可能で，引き続き直行する方向（z軸方向）に単純せん断変形を加えることができる．本報では，上述のように，解析モデルの形状をLx≒Lzに選んでおり，工学ひずみγ=2を加えるのに必要な上下（左右）のレプリカセルの変位量はそれぞれLx（≒Lz）になる．

Fig.3. Lees‑Edwards periodic boundary condition
2.5. 加工ルート
本研究では，準3次元条件を用いているため，押し込み加工ごとの被加工材の可能な回転は，Y軸まわりの90°回転，及び，Z軸まわりの180°回転の組み合わせのみである．なお，X軸まわりの180°回転は他の回転と一致する．そこで，ここでは表2に示すルートaからルートdの4つの加工ルートを考える．表中では，加工ルート名の後の括弧内に，ECAPの加工ルートを表すのに用いられる表記法13)を参考にしたものを併記している．つまり，記号 C はプロセスごとに被加工材を180°回転し， BC は90°回転し， BA は90°回転と‑90°回転を交互に加え，記号 A では被加工材を回転しないことを意味する．またγzxとγxzは被加工材を基準とした座標系o‑xyz における単純せん断ひずみを表し， + が工学ひずみγ=2，‑がγ=‑2を意味する．ルートa(A)とb(θZC)はECAPにおける代表的な加工ルートであり，ルートc(θYBC)とd(θZC‑θYBC)は，Y軸まわりの回転が必要であり実用的な加工ルートではないが，前述の理由から比較のために解析する．本報ではn回目の押し込み加工を第nプロセスと呼ぶこととし，ここでは4プロセスの加工（全工学ひずみ8）を加える．第3プ
ロセス以後は表中に示した第1，第2プロセスと同じ変形の繰り返しになる．各プロセス間では5psの緩和計算を行う．なお，せん断ひずみ速度は γ&= 2 × 1010s‑1とし，せん断変形の途中とプロセス間の緩和計算では，計算セルの寸法を一定に保つ．また，温度は300Kの一定になるように制御する．

３． 解析結果と考察
3.1. せん断応力の変化
図4に変形過程のせん断応力τzxと全付加ひずみγtの関係を示す．ここで，応力は個々の原子に対して有限温度での原子応力19)を評価し，それを全原子で平均化することで求めている．
図中では，各押し込み加工後の緩和計算（5ps）中の応力変化も示しており，緩和計算後の状態を，全付加ひずみに（ '）を付加して示している．また，今後，各プロセスのひずみをγによって表すことにする．前述のように全ての加工ルートで，第1プロセスは全く同じであり，第2プロセス以降に差異が現れる．なお，ルートaとbでは，変形条件から明らかなように，せん断応力の符号がプロセスごとに逆転する．図より，第1プロセスを除いて応力変化はほぼ周期的になっているのがわかる．それぞれの変形ルート間の違いを比較すると，
ルートcとdでは，各プロセスの最初に大きな応力を発生し，その後，ほぼ一定になることがわかる．一方，第3プロセスでの応力変化から比較できるように，ルートaでは，加工プロセスの最初に一度，応力のピークが現れ，γ≒0.5に，さらに大きいピークが見られる．
その逆に，ルートbでは，各プロセスの最初にピークが現れた後，ひずみγ≒0.5にかけて，大きく減少するのがわかる．また，各プロセスの最初に見られる応力のピークの大きさを比較すると，ルートaとbではルートcとdに比べて小さくなっている．これは，後に示すようにルートaとbでは，各加工プロセスの最初に容易に双晶変形を生じることが原因である．

3.2. 欠陥原子の割合の変化
図5にbcc構造以外の構造を持つ原子（欠陥原子と呼ぶ）の全原子に占める割合fdと全付加ひずみγtの関係を示す．ここで，原子構造は，CNA（Common Neighbor Analysis）14)を用いて評価している．ルートbでは，プロセスごとに欠陥原子の割合が徐々に減少し，それ以外のルートでは，第1プロセスを除いてほぼ周期的になる．ルートaとdでは，ひずみ0.5から1.2の間で欠陥原子の割合が最も高くなり，その後，減少する．その逆に，ルートbでは，欠陥原子の割合がγ≒0.8で最も少なくなり，その後，少し増加する．ルートcでは，他のルートと比べて変化が小さい．4プロセス後の最終状態を比較すると，ルートaにおいて最も欠陥原子の割合が高く，効率的に結晶内に欠陥が導入されているのがわかる．ルートbの欠陥割合が最も低くなり，ルートcとdではほぼ同じで中間的である．加工ルートによるせん断応力と欠陥原子の割合の変化が異なる原因については，次節で内部構造変化と関連させて明らかにする．

3.3. 内部構造変化
図6に第1プロセスにおける内部構造変化を示す．また，ルートa，b，c，dの第2プロセス以後の内部組織変化をそれぞれ図7，図8，図9，図10に示す．図中では，欠陥原子を黒色で表示し，bcc構造を有する原子に対しては，被加工材料の軸方向（図1(a)中のX軸方向に対応する方向）の結晶方位を基準としてRGBカラーを用いて色付けして表示している．この際，<0,0,1>，<0,1,1>，<1,1,1>が材料軸の方向を向いている結晶の色をそれぞれ赤，緑，青色としている15)．なお，ルートaとbでは，プロセスごとに材料の回転により結晶方位判定の軸が90°変化するため，図6(e)から図7(a)のように色が大きく変化する．図6より，変形を加えてすぐに多くの結晶粒内で双晶変形が発生し（図6(a)），初期結晶方位はすぐに消失する．個々の粒内で一通り双晶変形が終わると，その後は転位運動と粒界すべりによる変形が支配的となる．これは，せん断方向に双晶変形を一度生じると，同じ方向が逆双晶方向になり，双晶変形を非常に生じ難くなるためである．なお，内部組織が最も微細なのはγ=1.04（図6(b)）からγt=1.36（図6(c)）であり，図5において欠陥原子の割合が最も高くtなっている範囲に一致する．最終的には，材料軸の方向に<0,1,1>から<1,1,1>を持つ原子（図中で緑色から青色で表示された原子）が最も多くなる．図6(d)の変形を加えた直後の状態と図6(e)の緩和計算後の状態を比較すると，緩和により粒内の転位が減少し，非常に小さい粒が消滅するのがわかる．次にそれぞれのルートでの第2プロセス以後の内部構造変化を見ていく．ルートaではせん断変形の方向が変化すると，すぐに図7(b)のように双晶変形が現れる．最初に現れる双晶変形の方向はx 軸方向に近く，せん断変形の方向とは異なる．これはせん断ひずみと共役な応力成分が第1プロセスとは逆の方向，つまり順双晶方向に加わることが原因であると考えられる．さらに変形が加わると，図7(c)のようにz 軸方向に双晶変形が発生し，その前に発生した双晶境界を横切るため，図7(d)のように格子状の内部構造が現れる．この双晶変形の交差が，図4に示したせん断応力τzx−全付加ひずみ
γt関係において，ルートaのみ各プロセスで応力のピークが二度現れる原因である．ルート
aの第2プロセスでも，γ≒1（図7(d)）で様々な方位を持つ結晶が含まれており，最も組織が微細である．さらに変形が加わると，第1プロセスと同じように全体的に材料軸の方向に<0,1,1>から<1,1,1>を持つ原子が増加する．第3プロセス以後は同様の変形の繰り返しになるが，各加工プロセス後の組織（図7(h)，(i)，(j)）を比較すると，付加ひずみが大きくなる程，結晶粒子数が増加するのがわかる．ルートbではせん断変形が逆になるとすぐに図8(a)のようにx軸と並行な双晶変形が現れる．さらに変形を加えていくと，図8(b)，(c)のように第1プロセス後半と同様の組織になる．図8(d)，(e)からわかるように，ルートbでは全付加ひずみが大きいほど，結晶粒子数が減少しており，このような繰り返し変形は微細化プロセスとして不適切であることがわかる．なお，より長時間の緩和計算を実施すると，図8(e)中に見られる微小な粒も消滅すると考えられる．ルートcでは常に加わるせん断変形の方向が同じであるため，細かい結晶粒の結合，転位運動による粒界の形成が多少見られるものの，他のルートほど顕著な構造変化は見られない．初期結晶粒が小さいこともあり，すぐに定常状態に達するのがわかる．ルートdでは，第2プロセスのせん断変形の方向がルートbと逆であるため，せん断ひずみと共役な応力成分が第1プロセスにおける単純せん断と同じ方向（逆双晶変形の方向）に加わる．そのため，図7(b)のようにx 軸方向の双晶変形はほとんど生じず，z軸方向の双晶変形が支配的である．4つの加工ルートにおける最終的な組織（図7(j)，図8(e)，図9(e)，図10(e)）を比較すると，ルートaが最も組織が微細であり，双晶変形が交差するメカニズムが非常に有効であることがわかる．また，Y軸まわりの材料回転が必要なルートcとdは，結晶粒微細化効率においても劣る．なお，強変形中にまれに微小なボイド状の欠陥が生成するのが確認されたが，それらはすぐに消滅し，成長することはなかった．
４． まとめ本研究では，L‑E周期境界条件を用いたMDシミュレーションにより，金型角度がΦ=90°の場合のECAPと等価な強変形を実現し，準3次元条件下で実現可能な4通りの加工ルートにおける内部構造変化の解析を実施した．その結果，加工プロセスごとに材料に回転を加えない最も単純なルートa(A)において，最も効率的に組織が微細化されることがわかった．これは，単純せん断と共役な応力成分によって生じる双晶変形と，せん断変形の方向の双晶変形が交差するためである．またルートaの各加工プロセスで最も組織が微細になるのは，γ≒1の時であることがわかった．準3次元条件下の計算ではECAPにおける他の代表的加工ルートであるルートθZBCとθZBAによる内部構造変化の解析を行なうことができない．より大規模な計算が必要になるが板厚方向の寸法を十分に厚くし，L‑E周期境界条件と周期セルの変形を組み合わせることで，それらのルートにおける組織変化を解析することも可能であり，今後の課題としたい．謝辞
本研究に対しご援助頂いた財団法人 JFE21 世紀財団に深甚なる謝意を表します。
参考文献
1. E. O. Hall, Proceedings of the Physical Society, Section B, 64, pp.747‑752 (1951).
2. N. J. Petch, Journal of the Iron and Steel Institute, 174, pp.25‑28 (1953).
3. A. Makino, A. Inoue and T. Masumoto, Material Transactions JIM, 36‑7, pp.924‑938(1995).
4. C. Suryanarayana and F. H. Froes, Journal of Material Research, 5‑9, pp.1880‑1886(1990).
5. R. Z. Valiev, A. V. Korznikov and R. R. Mulyukov, Materials Science and Engineering

A, 168‑2, pp.141‑148 (1993).
6. M. Nemoto, Z. Horita, M. Furukawa and T. G. Langdon, Metals and Materials

International, 4, pp.1181‑1190 (1998).
7. Z. Horita, D. J. Smith, M. Furukawa, M. Nemoto, R. Z. Valiev and T. G. Langdon,

Journal of Material Research, 11, pp.1880‑1890(1996).
8. Y. Saito, H. Utsunomiya, N. Tsuji and T. Sakai, Acta Materialia, 47, pp.579‑583(1999).
9. 北川浩，北村隆之，渋谷陽二，中谷彰宏,

初心者のための分子動力学法, (1997)

養賢堂.
10. 上田顯,

分子シミュレーション

11. 佐藤明,How to 分子シミュレーション, (2003) 裳華房.
, pp.25‑27 (2004) 共立出版.

12. G. Simonelli, R. Pasianot and E. J. Savino, Materials Research Society Symposium

Proceedings, 291, pp.567‑572 (1993).
13. M. Furukawa, Z. Horita and T. G. Langdon, Materials Science and Engineering, A332,pp.97‑107 (2002).
14. J. Dana Honeycutt and Hans C. Andersen., The Journal of Phisical Chemistry, 91‑19,pp.4950‑4963(1988).
15. 松本龍介，北川浩，中谷彰宏, 材料, 52‑3, pp.235‑240 (2003).

本研究に関連する発表
論文
1. 松本龍介・林田登志男・中垣通彦， 準３次元条件下での ECAP によるα鉄における結晶粒微細化過程の分子動力学解析 ，材料．（投稿中）
国際会議での発表

1. R. Matsumoto, T. Hayashida & M. Nakagaki, “Molecular Dynamics Analysis on
Initial Texture and Processing Route Influence on Grain Refinement Behavior of α-Fe by Equal Channel Angular Pressing”, The 8th Asia-Pacific Symposium on
Engineering and its Applications (AEPA 2006), Sep. 2006. (will be presented)．
国内学会での発表
1. 林田登志男・松本龍介・中垣道彦， 強せん断変形による多結晶金属の内部構造変化の分子動力学シミュレーション ，日本機械学会（M&M 材料力学カンファレンス），No.05-9，(2005-11)，pp.559-560．
2. 林田登志男・松本龍介・中垣道彦， ECAP による組織変化の分子動力学シミュレーシ
ョン ，日本機械学会（第１８回計算力学講演会），No.05-2，(2005-11)，pp.467-468．
3. 林田登志男・松本龍介・中垣道彦， ECAP によるα鉄の結晶粒微細化挙動に与える初期組織と加工ルートの影響の分子動力学法による検討 ，日本機械学会（九州支部第５９期総会講演会），(2006-3)．（講演予定）