塑性加工による集合組織形成に関する分子動力学解析
山形大学 大学院理工学研究科 機械システム工学分野
教授 上原 拓也
（平成 26 年度一般研究開発助成 AF-2014018）
キーワード：微視組織，分子動力学法，大変形

１．研究の目的と背景
圧延などの大きな塑性変形を伴う加工中には集合組織を用いて数値積分することによって，すべての原子の運動を解く．

が形成され，このような集合組織を含む微視組織が材料の加工性や力学特性などに大きく影響することが知られている．そのため，逆に微視組織を制御することによって所d ri  1 Fi , Fi   f ij   ij
mij ij i drij(1)

望の特性を得ることが可能であり，形成される微視組織を制御する様々な技術が開発されつつある．また，加工プロ
セスや機械構造物の変形・強度予測には，計算機シミュレーションが不可欠となっており，計算精度の向上のためには，詳細な材料データが必要である．このような材料データが微視組織に依存し，さらに材料組成や温度，応力，ひずみ速度など，様々な条件にも依存するため，これらのデータを収集するには，膨大な実験が必要となる．微視組織の影響を考慮した力学解析としては，結晶塑性論などを用いて加工中の変形と微視組織の関係を考慮した数値解析手法の構築が進められているが，これらを用いた計算精度の向上に必要なだけの十分なデータを実験によって収集するのは困難であるのが現状である．そのため，材料データそのものを計算機シミュレーションで求めようとする研究が進められている．例えばフェーズフィールドモデルを用いると，凝固組織や多結晶体などの微視組織が再現でき，そこから材料特性の微視組織依存性を推測することが可能である．しかしながら，力学特性に大きく影響する結晶方位変化や粒界挙動は，本質的には結晶構造や原子配列の変化を伴う現象であり，これらを明らかにするためには，原子スケールの解析が不可欠である．
そこで本研究では，塑性加工中に生じる集合組織の形成過程を原子スケールのシミュレーションによって再現し，集合組織形成のメカニズムを明らかにすることを目的と

ここで，ri は原子 i の位置ベクトル，mi は原子 i の質量，
Fi は原子 i に働く力であり，これは周囲の原子 j から働く原子間力 fij の和で表されるとする．また，fij は二体間ポテンシャル関数ij で表されるとものとし，このポテンシャル関数としては，以下で表される Lennard-Jones 型を用いる．

ij = 4D {(R/rij)12-(R/rij)6}(2)

ここで，D と R は材料パラメータであるが，無次元化した解析を行うことによって，これらの数値は定性的には影響しない．
2・2 解析モデル研究で用いた計算モデルを図１に示す．計算モデルは３次元であるが，図のような２次元的な領域パターンを作成する．すなわち，x-y 平面において正方形の計算領域を４つの領域に分割し，互い違いになるように，２列ずつ配置する．これらの領域が結晶粒に相当し，それぞれに I～IVの識別記号を付す．このとき，負荷方向（y 軸方向）に対して I と II，および III と IV が直角方向に並ぶモデル（Model A）と，平行に並ぶモデル（Model B）を作成する．今回用いるポテンシャル関数では，結晶構造として fcc 構造が安定であることから，初期配置として，すべての原子を fcc のする．将来的には，既に存在する集合組織が加工特性に与える影響についても原子スケールのシミュレーションから明らかにすることが可能とすることを目的とするが，本Compression

Compression
III

研究では，その基本的な計算モデルの提案とその妥当性のIII検討を行う．さらに将来的には，所望の加工性や強度特性IVIIIII
IVをもつ材料を得るための微視組織形態を計算力学モデルから推定し，その組織を得るための加工条件を逆解析的に求める手法を構築することを目指した研究を推進する．

２．解析方法とシミュレーションモデル
2・1 解析手法
本研究では，古典的分子動力学法を利用する．基礎式は以下に示す Newton の運動方程式であり，これを Verlet 法

Grain III

Grain Iyz

IIIx(a) Model A
図１．計算モデル(b) Model B

格子点に置いた後，図の領域ごとに異なる方位となるように回転して初期配置とする．境界条件は周期境界条件とし，図の奥行き方向（z 軸方向）にも周期境界条件を与える．モデルサイズは，x, y, z 方向にそれぞれ 50, 50, 5 単位セルとする．初期配置の作成において，領域の境界，すなわち結晶粒界付近では，極端に近すぎる原子ペアが(a) 12000 step(b) 20000 step(c) 24000 step(d) 26000 step(e) 30000 step(f) 36000 step(g) 42000 step(h) 48000 step(i) 54000 step(j) 60000 step(k) 64000 step(l) 70000 step

生じた場合には，片方の原子を削除することで，大きな反発力によって計算が不安定となることを避けている．
この操作のため，原子数はモデルによってわずかに異なるが，いずれのモデルでもおよそ 50000 である．
2・3 解析条件
負荷条件としては，15000 ステップの間，無負荷状態で緩和計算をした後，y 軸方向の長さを一定速度で変化させる．本報では，Ly = 25.0 だけの圧縮を 50000 ステップの間に負荷する．緩和状態での y 軸方向の長さが Ly = 73.0程度であるので，およそy = 0.34 の圧縮ひずみに相当する．
この間，x, z 方向には，直方体形状を保ったまま，辺の長さを変化させることによって，応力の各方向の垂直成分が０となるように制御する．圧縮過程の終了後，序か過程として，5000 ステップの間，y 軸方向の応力成分が０となるように制御した計算を行う．

図２．原子配列の変化（モデル A1）

図１に示した結晶粒 I～IV の初期方位としては，表１に示す３通りについて計算した結果を示す．ここで，26.6°

ら，図 2 (c) に示す 24000 ステップでは，右上と右下の結は tan 1/2，18.4°は tan 1/3 に相当し，規則的な対応粒晶粒界付近に，濃い青色で示される領域が形成される．こ界が得られる方位として設定した．の濃淡は，初期配置では [001] 方向を z 軸にとっていたた-1-1

め，奥行き方向の原子がすべて重なって投影されていたの表１．初期結晶方位に対し，この領域ではずれが生じたことを表している．つまり，この領域で結晶方位の変化がおきたことを示していID
A1/B1A2/B2A3/B3Grain NoIIIIII
IV+26.6° -26.6° +18.4° -18.4°+26.6° -18.4° +18.4° -26.6° Illustration of the+26.6° +18.4° -26.6° -18.4° case for 26.6°

る．この方位変化の領域は速やかに結晶粒内に広がるとともに，他の結晶粒でも方位変化が生じ，図 2 (d) に示す
26000 ステップではすべての結晶粒で方位変化が見られる．この方位変化は，図 2(c)で見られたように，結晶粒界，とくに粒界３重点が起点となることから，同じ結晶粒でも異なる起点からの領域が発生し，それらが衝突する面が形３．計算結果と考察成される．その結果，もともとの結晶粒内に，複数の結晶3・1 圧縮変形方位領域が存在し，いわゆるサブグレイン領域が形成されまず，モデル A1 での計算結果を示す．図２は，圧縮変たとみなすことができる．ただし，このとき形成される新形過程における変形と原子配列の変化を表している．個々

しい粒界は，初期に設定した大計画粒界とは異なり，安定の原子のポテンシャルエネルギーによって色分けをしてではなく，一時的に形成しただけの場合もあり，図 2 (e)表示しており，青が低く，緑，黄，赤となるにつれて高いに示す 30000 ステップでは，左下の結晶粒（Grain I）では，エネルギーを表している．結晶粒内では，原子は fcc 構造

再び龍内が均一な方位にそろっている．それに対し，右下の格子点に規則的に配置しており，安定な状態であるのにの結晶粒（Grain II）では，２つのサブグレインが残留し対し，結晶粒界では原子配置が不規則となり，エネルギーている．

が高まるため，この図では，青以外の原子を見ることで結さらに圧縮変形を加えていくと，図 2 (f) に示す 36000晶粒界が表現されている．図 2 (a) は，圧縮負荷を加えるステップでは，もともとは縦方向（y 軸方向）であった結前の緩和状態であり，初期配置をして与えた４つの結晶粒晶粒界が，押しつぶされるように斜め方向に傾き，結晶粒が，ほぼ正方形形状を保っていることがわかる．このモデの形状が大きく変化することがわかる．その後は，変形がルに対し，y 軸（図の上下方向）に圧縮負荷を加えると，進むにつれ，結晶粒内での方位変化および新たな結晶粒界

y 軸方向の長さ Ly が短くなるが，20000 ステップでは，原の形成が断続的に発生し，細かい結晶粒も観察されるよう子配列に大きな変化は見られない（図 2 (b)）．しかしながになり，図 2 (j) に示す 60000 ステップになると，もはや初期の結晶粒はほとんど判別できなくなっている．なお，の方位変化が見られるものの，大きな変化はないことから，圧縮変形が進み，y 軸方向の長さがかなり短くなるこのあ

降伏による原子配列の乱れが完了した後，原子配列は一旦たりから，y 軸方向で隣接する（上下に並んでいる）結晶安定化し，再び弾性的に変形しているとみなすことができ粒の結合がみられる．これは，y 軸方向にも周期境界条件る．

をかけているため，周期性を満たす方位として，制限が強くかかることが影響している．

その後，30000 ステップを過ぎたあたりで，応力は再びピークを示し，急激な応力低下が見られるが，このときの

図 2 (k)に示す 64000 ステップの後，65000 ステップから応力値は，一度目のピーク値より小さな値である．原子配は除荷し，y 軸方向の圧縮応力を０にする．その結果，弾列では，図 2 (f) のように，ピーク前から大きな変化を示性ひずみ分の回復が起こるが，微視組織については，除荷していることが確認できる．その後，応力はやや上昇傾向によってはほとんど変化せず，最終状態の 70000 ステップは見られるものの，最初の降伏前や，２度目のピーク前に（図 2(l)）で計算を終了する．

見られたような弾性的な応力上昇はみられず，明確な応力3・2 応力変化ピークや急激な応力低下も見られなくなる．これは，図 2

図２に示した過程における負荷方向の応力成分 y の変(f)～(k)で確認できるとおり，この段階では結晶方位の変化を表したのが図３である．ただし，圧縮を正としており，化や，新たな粒界の形成が断続的にあちこちで発生してい図中には y 軸方向の辺の長さも示している．横軸は時間スることに起因する．なお，65000 ステップで除荷する際にテップであるが，辺の長さの変化率は一定であるため，圧は，ほぼ瞬間的に弾性分のひずみが回復して，応力が０に縮負荷中となる 15000～65000 ステップにおいては，横軸なっていることがわかる．

を公称ひずみと読み替えることができ，定性的には応力ひ3・3 結晶方位の影響ずみ線図と等価である．

次に，結晶粒の配置はそのままにして，各結晶粒の方位圧縮負荷を開始した 15000 ステップから 24000 ステップを変えたモデル（Model A2）での計算結果を示す．図４は付近までは，応力は単調に増加する．傾きはやや鈍る傾向

原子配置の変化を表しており，応力負荷前の緩和計算では，にあるが，ほぼ直線的な増加とみなすことができる．前節

モデル A1 のときと同様に，正方形状の結晶粒が保持されで説明したとおり，この期間には原子配置や結晶粒形状にている（図 4 (a)）．この状態に圧縮負荷を加えていくと，はほとんど変化が見られなかったことも踏まえ，この段階はじめは弾性変形を生じるのみで，結晶粒や原子配置に変が弾性変形領域であるとみなされる．その後，応力は化は見られない（図 4 (b)）が，24000 ステップになると24000 ステップあたりで急激に減少する．図 2 (c), (d) に１つの結晶粒（Grain III）内に，方位変化を起こした領域示したとおり，ここで原子配列が大きく乱れており，非弾が形成される．この変化の起点となったのは右下の粒界３性変形が起こったといえる．マクロな応力ひずみ線図でい重点であり，ただちに結晶粒内の全領域に広がっていく

う降伏点や降伏応力との対応についての詳細な検討は避けるが，本報ではこの点を「降伏点」とよぶこととする．
降伏によって急激に低下した応力は，再び上昇をはじめる．このときの上昇率は，降伏前と同じ程度であり，弾性変形的な特徴を示している．図２ (d), (e)に示した
26000 および 30000 ステップの原子配列では，結晶粒内でLy

Stress in y y(b) 20000 step(c) 24000 step(d) 26000 step(e) 30000 step(f) 36000 step(g) 42000 step(h) 48000 step(i) 54000 step(j) 60000 step(k) 64000 step(l) 70000 step7065260
155050Length in y LyPy3(a) 12000 step15000 25000 35000 45000 55000 65000
Time step k

図３．応力変化（Model A1）

図４．原子配列の変化（Model A2）

(a) 12000 step(b) 20000 step(c) 26000 step(d) 28000 step(e) 30000 step(f) 36000 step(a) A1(b) A2(c) A3(d) B1(e) B2(f) B3

図６．結晶粒分布の比較（48000 step）

4(g) 42000 step(h) 48000 step

A1A2A3B1B2B3(i) 54000 step(j) 60000 step(k) 64000 step(l) 70000 step

図５．原子配列の変化（モデル B1）（図 4 (d)）．この状態では，他の結晶粒内にも，新たな結晶粒界が形成され，サブグレインとよべるような領域が見Stress in y y32100 0.05 0.1 0.15 0.2 0.25 0.3
Strain in y yられる．このような一連の大きな変化は，一旦，収束し，26000 ステップから 30000 ステップの間にはあまり変化が図７．応力ひずみ線図の比較見られない．その後，36000 ステップまでには再び結晶粒の分布が大きく変化し，微細な結晶粒も形成されていく．

お，60000 ステップ以降を表す図 5 (j)～(l) において，結このように，定性的にみれば，結晶方位の相違にはよらず，晶粒ごとに濃淡が違って見えるのは，3.1 節で説明したと前節に示した Model A1 の結果と同様な挙動を示すことがおり，z 軸方向への原子の重なりの影響であり，ここでは確認できる．

偶然的に方位がそろったためである．

3・4 結晶粒配置の影響

3・5 考察次に，3.1, 3.2 節に示したモデル（Model A1）に対し，最後に，表１に示した６モデルのすべてについて，計算

各結晶粒の結晶方位はそのままにして，結晶粒の配置を変結果をまとめて示す．図６は，48000 ステップにおける原えたモデル（Model B1）に対する計算結果を図５示す．応子配列であり，図７は圧縮および除荷過程における応力ひ力を付加しない状態（図 5 (a)）では，この場合も正方形ずみ線図である．

状の結晶粒が保たれており，20000 ステップあたりの弾性原子配列の変化については，どのモデルにおいても，初変形範囲では，結晶粒分布に変化は見られない（図 5 (b)）期の結晶粒界に加え，多くの粒界が形成されることによっが，図 5 (c) に示す 26000 ステップになると，結晶方位がて，結晶粒の微細化が起こっている．結晶粒の形状や結晶変化した領域が現れ，28000 ステップにかけて，結晶粒全方位に関して，負荷方向との関連性や，新たな結晶粒方位体に広がっていく（図 5 (d)）．ここで，図 3 に比べて特徴とオリジナルの結晶方位の関係などについて，詳細な検討的な点は，負荷方向（y 軸）に対して垂直な方向の粒界はをする必要がある．塑性加工などの大変形に伴う集合組織そのままの方向を保持しているのに対し，負荷方向に平行としては，このような関係を明らかにする必要があり，今な粒界が斜めに傾く点である．また，その傾く角度が，互後の継続研究課題としたい．

い違いになるため，正方形であった結晶粒が，六角形状にまた，応力ひずみ線図においても，すべての場合におい変形している点である．また，結晶粒内に新たに形成するて，初期には応力がほぼ線形に上昇した後，応力が急激にサブグレインの粒界は，総じて，縦方向に現れる傾向が見低下し，再び応力が線形的に上昇した後，２度目の応力ピられる．

ークを迎えて応力が急激に低下，その後は断続的に応力のその後，圧縮ひずみが大きくなるにつれ，微細な結晶粒上昇と低下が見られるという定性的な傾向が共通して見が現れていく傾向は，これまでのモデルと同様である．なられた．ただし，降伏応力は，最大値は 3.2 程度，最小値は 2.8 程度である．また，２度目の応力低下以降の挙動に

参考文献

ついては，図３に示した Model A1 とは異なり，Model A3,

1) T. Uehara: Molecular dynamics simulation of stick-slip

B3 などでは，明確な応力ピークと急激な応力低下を周期

friction on a metal surface, Appl. Mech. Mater., 459 (2014),的に繰り返し，Model B3 では４度目のピークが最大の応

pp. 26-33.

力を示すなど，特徴的な挙動を示している．このような特

2) T. Uehara: Molecular dynamics simulation of tensile

徴が現れる要因についても，結晶方位や粒界と負荷方向のproperties of nano-layered materials, Adv. Mater. Res., 741

関係などと関連づけて考察することによって，今後，明ら(2013), pp. 79-83.

かにしていく予定である．

3) T. Uehara and H. Suzuki: Numerical simulation ofhomogeneous

４．結言polycrystalline

grain

formationusing

multi-phase-field model, Appl. Mech. Mater., 197 (2012), pp.

本研究では，分子動力学法を用いて，圧縮の大変形を受

2610-2614.

けるときの微視組織変化をシミュレートし，初期の多結晶

4) T. Uehara: An atomistic study on the slip deformation

体における結晶粒内に結晶方位変化が起き，サブグレイン

mechanism of crystalline materials using a weak-plane領域が形成されることを示した．本研究では，詳細な結晶

model, Appl. Mech. Mater., 197 (2012), pp. 2070-2075.

方位や，負荷方向との関連性などについては，十分な解析

5) T. Uehara: A phase field modelling for multi-scaleを行うまでには至らなかったが，今回提案したモデルの有deformation mechanics of polycrystalline metals, Procedia

効性を示すことができた．今後は，より詳細な解析を行う

Eng., 10 (2011), pp. 1779-1784.

とともに，変形負荷前後の力学特性変化を評価することに

6) T. Uehara: Numerical study on the evolution of stress

よって，塑性加工による微視組織変化と力学特性変化を統

distribution in cellular microstructures, Int. J. Mater. Res.,

一的に解析できる解析手法の構築に向け，基礎的なモデルを示すことができたといえる．

101 (2010), pp. 486-491.
7) T. Uehara, C. Asai, and N. Ohno: Molecular dynamics
simulation of shape-memory behaviour using a multi-grain

謝 辞
本研究は，公益財団法人

model, Modell. Simul. Mater. Sci. Eng., 17 (2009), #035011.

天田財団の一般研究開発助成

8) T. Uehara: Molecular dynamics simulation of the variation in（AF-2014018）の支援を受け実施した．ここに深く感謝

the microstructure of a polycrystalline material under tensile

申し上げます．

load, Key Eng. Mater., 748 (2017), pp. 375-380.
9) T. Uehara: Change in microstructure of a polycrystalline
material under compression, Proc. XIV Int. Conf. on Comp.
Plasticity (COMPLAS XIV), (2017), pp. 106-113．