Multi-Phase-Field 法の並列 GPU 計算によるフェライト変態の高精度予測研究代表者

東京農工大学大学院工学研究院 先端機械システム部門准教授 山中晃徳

1. 緒言
高強度鋼の機械的特性は、材料に内在するミクロ組織の形態や分布状態に大きく左右される。たとえば、代表的な高強度鋼である Dual-Phase(DP)鋼においては、第二相のマルテンサイト相の体積分率や結晶粒径が同じであっても、マルテンサイト相の分布が異なれば、強度や延性が異なることは良く知られている。そのため、高強度鋼のより高精度な材質制御を目的として、ミクロ組織の形態や分散状態と機械的特性の関係を明らかにしようとする研究が盛んに行われている。また最近では、ミクロ組織の 3 次元観察および定量評価が世界的に注目されている(1)~(3)。足立ら(4)~(6)は、機械的研磨によるシリア
ルセクショニングと光学顕微鏡または SEM-EBSD を用いて鉄鋼材料中の相変態組織の3 次元像を取得し、その組織形態を位相幾何学や微分幾何学によって数学的に定量評価（計量形態評価）することを可能としている。さらに彼らは、DP 鋼のフェライト-マルテンサイト二相組織の 3 次元形態をガウス曲率などの数学的な形態特徴値により定量評価し、それらの組織特徴値と DP 鋼の機械的特性との関係を解明しようとしている。
このような実験的な手法は非常に有効であるが、ミクロ組織の三次元形態の評価や機械的特性との関係をすべて明らかにすることは、当然ながら膨大な時間と金銭的コストを要する。特に、熱処理におけるミクロ組織形態の時間変化やその変化が機械的特性に及ぼす影響を明らかにすることは、実用的にも学術的にも極めて重要であるが、実験的な手法のみではほぼ不可能であろう。
本研究では、基礎的ではあるが高強度鋼の複相組織形態を理解するうえでも必要不可
欠な Fe-C 合金で生じるオーステナイト→フェライト（γ→α）変態挙動の高精度予測およびフェライト相形態変化の定量評価を実現するための数値シミュレーション法を確立することを目的する。本報告では、本研究で確立したマルチフェーズフィールド(Multi-Phase-Field; MPF)法による Fe-C 合金のγ→α変態シミュレーションとシミュレーション高速化のための Graphics Processing Unit (GPU)による並列計算、さらに MPF シミュレーションで得られるフェライト相の計量形態評価法について報告する。
2. シミュレーション方法
2.1 Multi-Phase-Field 法
本研究では、Fe-C 合金で生じるγ→α変態挙動の数値シミュレーション法として、Steinbach ら(7)の MPF 法を用いる。MPF 法は、ミクロ組織の全自由エネルギーの最小化原理に基づき、ミクロ組織の形成過程を数値解析する方法であり、モンテカルロ法や－ 79 －

セルオートマトン法の問題点を克服した、非常に強力なシミュレーション法として注目を集めている(8)~(9)。以下で、本研究で用いた MPF モデルを説明する
2.1.1 秩序変数と全自由エネルギー
MPF 法では、対象とする現象を表現するために必要十分な物理量を座標 x、時間 tの関数である「場の秩序変数」として定義する。Fe-C 合金における γ→α 変態では、フェライト相の成長と炭素原子の拡散が生

Fig.1 Definition of i(x,t)

じるため、各相(結晶粒)の存在確率を表す

フェーズフィールド変数i(x,t)と炭素濃度変数 C(x,t)を秩序変数として定義する。フェーズフィールド変数i(x,t)は、結晶粒 i の存在確率を表すため、Fig.1 に示すように結晶粒i 内ではi(x,t)=1、他の結晶粒内ではi(x,t)=0 の値となり、結晶粒 i と他の結晶粒との界面領域で 0 から 1 に滑らかに変化するものと定義される。また、フェーズフィールド変数は任意の座標において、次の条件を満足する必要がある。
N x,t  1i1(1)i一方、炭素濃度 C(x,t)は、結晶粒 i における局所的な炭素濃度(以後、これを局所濃度と呼ぶ)ci(x,t)を定義し、各結晶粒のフェーズフィールド変数の重み付け和として、次式で定義する。
Cx, t  N x, t c x, t i 1i(2)i

なお、後述するように、局所濃度は Fe-C2 元系平衡状態図を線形近似することにより求める。

MPF 法では、「ミクロ組織の形成は、系（ミクロ組織）の全自由エネルギーが単調減少するように進行する」との仮定に基づいている。本研究では、系の全自由エネルギーを、次式のギンズブルグ-ランダウ型のギブスの自由エネルギー汎関数として定義する。
G N N  aij
 N N
 N       i x, t   1dVW


g



ijij i jchemV i 1 j i 1 2
 i 1 i 1 j i 1  (3)ここで、右辺の被積分項の第一項は、勾配エネルギーと呼ばれ、ミクロ組織中に界面が存在することで、存在しない場合に比べて過剰に生じるエネルギー量を表す。第二項は、ダブルウェルポテンシャルと呼び、異なる結晶粒が同じ座標に存在することを防ぐペナルティー項である。また第三項は、化学的自由エネルギー項である。第四項は、式(1)の条件を満たすための項であり、はラグランジュの未定数定数である。さらに、式－ 80 －

の aij と Wij は、それぞれ勾配係数、ポテンシャル障壁高さと呼ばれるパラメータであり、2.1.3 で説明するように、物性値である界面エネルギーや界面幅に関係づけられる。
2.1.2 秩序変数の時間発展方程式
MPF 法においては、フェーズフィールド変数i(x,t)の時間発展方程式は、次式で表される。
 G G 
 i2 N
M ij  i  j t
N j i(4)

ここで、 M ij はi(x,t)の増減を決めるフェーズフィールドモビリティーと呼ばれ、一般には界面の易動度と関係づけられる。また、式(4)に式(3)を代入することで、ラグランジュの未定乗数は相殺され、次式が得られる。
N
g
G G
1
  g2 aik  a jk   k  Wik  W jk  k    chem  chem  j i  j k 1  2
  i(5)

さらに、式(5)を式(4)に代入すれば、次式に示すi(x,t)の時間発展方程式が導出される。
N N 1
 i
g2
  g
M ij   a ik  a jk  2 k  Wik  W jk  k    chem  chem  2N jit ij  k 1  (6)

ここで、化学的自由エネルギー密度 gchem が、フェーズフィールド変数、炭素濃度、温度の関数として与えられる場合、式(6)右辺第二項は具体的に計算できるが、本研究では次式で近似する。
gchem gchem
8i j Eiji j(7)

ここで、Eij は相変態の化学的駆動力の大きさを表す。
さらに、式(7)の右辺においては、N 個の結晶粒について、結晶粒間の相互作用(Σの計算)を計算しているが、N 個すべての結晶粒間の相互作用を計算するのは、計算時間の増大にもつながる。つまり、界面を構成している結晶粒についてのみ、相互作用の計算をすればよい。
そこで、各座標において界面を構成している結晶粒の個数を次式で計算する。
Nnx, t   i(8)i 1

ここで、�� は� � �� ≤ 1を満たす領域で 1、それ以外では 0 の値となるステップ関数である。したがって、 2 つの結晶粒が接する界面では n = 2、粒界 3 重点では n = 3 となる。これより、式(6)の N を n に置き換えることで、i(x,t)の時間発展方程式は、最終的に次式で与えられる。
n n 1
 i2
 8
M ij   a ik  a jk  2 k  Wik  W jk  k   i j E ij tn j i  k 1  2

 

(9)

一方、炭素濃度 C の時間発展方程式は、結晶粒 i からの拡散流束 Ji の重みづけ和として、次式で表される。

－ 81 －N NC    i J i      i DiC ci t i 1
 i 1(10)

ここで、‫ܦ‬௜஼ は結晶粒 i における炭素原子の拡散係数であり、結晶粒 i がα相かγ相かによって値が異なる。式(10)の説明を簡単にするため、以下では 2 つの結晶粒からなる系(N = 2)を考え、結晶粒 1 がα相、結晶粒 2 がγ相とする。このとき、式(10)は次式のように表すことができる。

C   1 D1C  c1   2 D 2C  c 2
t



(11)

ここで、結晶粒 1 と 2 がそれぞれα相とγ相であることを上付き添え字で表した。Fe-C 合金において、γ→α変態が生じると、α相で固溶できない炭素原子がγ相へと拡散し、濃度分配が生じる。このときのα相とγ相間の炭素濃度の分配は、次式の濃度分配係数 kij を用いて表現する。
C1 C eqkC 2
C eq

(12)

C 2 C eq
k 22      1
C 2 C eq

(13)

k12 

ここで、k12 は、γ相である結晶粒 2 を基準としたときのα相である結晶粒 1 の濃度比を表してఊఈ）の比率に等しいものとする。式と式を、式
おり、α相とγ相の平衡炭素濃度（‫ܥ‬௘௤ 、‫ܥ‬௘௤の炭素濃度の定義式に代入すれば、結晶粒 1(α 相)と結晶粒 2(γ 相)の局所濃度は次式で与えられる.

c1 

k12C
k121  k222

(14)

c2 

k22C
k121  k222

(15)



したがって、式(14)と式(15)を式(11)に代入すれば、炭素拡散方程式は次式で表される。


2
CC C    i Di k i 2  2t  i 1kj2j j 1(16)なお、各相における炭素原子の拡散係数は、次式のアレニウス型の式がよく用いられる。
 Q 
DiC  D0Ci exp i  RT 

(17)

2.1.3 パラメータと物性値の関係づけ
式(9)で示したi(x,t)の時間発展方程式の計算に必要なパラメータであるܽ௜௝ とܹ௜௝ は、定常状態における界面領域を考慮することにより、次式のように物性値と関係付けられる。

－ 82 －

4 ij

W ij aij (18)



8ij

(19)

2

ここで、ij は界面エネルギー、ߜは界面幅である。
さらに、式(9)のフェーズフィールドモビリティーに関しては、α/γ界面の易動度と次式のように関係付けられる。

2M ij8

M ij (20)

ここで、Mij が結晶粒 i と結晶粒ｊの界面の易動度であり、α/γ界面の場合のみ非零とし、アFig.2 Linearized phase diagram of Fe-C system

レニウス型の次式を用いる。

  QM M   M 0 exp RT (21)2.1.4 Fe-C2 元系平衡状態図との関係付け
式(9)の߶௜ の時間発展方程式の変態駆動力Eij と式(11)の炭素拡散方程式の濃度分配係数 kij を Fe-C2 元系平衡状態図を線形近似することで求める方法を説明する。

Fig.2 に、Fe-C2 元系平衡状態図の線形近似の概念図を示す。図中の実線は、CALPHAD法で得られる α 相と γ 相の化学的自由エネルギー関数に基づき共通接線則から求めた平衡状態図である。破線の直線は、この平衡状態図を線形近似したものである。このとき、ܶ௜ோ を平衡状態図を線形近似する際の基準となる温度とすれば、温度 T とα相およびγ相の平衡炭素௘௤௘௤

濃度‫ܥ‬ఈ 、‫ܥ‬ఊ の関係は、次式のように表される。
T  TR  m Ceq
R(22)
(23)

eq

T  T  m C

ここで、݉ఈఊ は線形近似されたα/α＋γ境界線の勾配、 ݉ఊఈ は線形近似されたγ/α+γ境界線の勾配である。この関係式より得られる、任意の温度 T における各相の平衡炭素濃度を用いることで、式(12)と式(13)より分配係数 kij を計算できる。

一方、変態駆動力Eij は、変態によって生じるエントロピー S と過冷度ΔTij の積で表す。
(24)
Eij  STij
ここで、過冷度ΔTij は、結晶粒 i と j がα相とγ相である場合のみ非零となり、次式で表される.
T T  T
2(25)

T

ここで、ܶఈఊ とܶఊఈ は式(22)と式(23)を用いて次式にて与えられる。
T  TR  m (ci  CR )(26)

T  TR  m (cj  CR )(27)－ 83 －

2.2 GPU による並列計算
γ→α変態により形成されるフェライト相の3次元形態評価を行うために、MPFシミュレーションの3次元計算を行わなければならない。しかしながら、MPF法は秩序変数と同数の複数の時間発展方程式の数値計算を行う必要があり、計算コストが大きい。そこで、計算機アクセラレータとして注目され

Fig.3 Domain decomposition for multiple-GPUているGPUを用いた計算の高速化をcomputing行う。さらに本研究では、現在のGPUのメモリ容量は最大で6GB程度(NVIDIA TESLA K20X)であることから、現実的な3次元空間におけるMPFシミュレーションを行うために、複数GPUを用いた並列計算法を確立した。
複数 GPU を用いて並列計算を行う際には、計算領域全体を複数の副領域に分割し、それぞれに 1 基の GPU を割り当てる。本研究では、Fig.3 のように計算領域全体を z 方向に 2 つの副領域(subdomain)に 2 次元分割し、2 基の GPU を用いて計算を行う。次に、各 GPU が担当する副領域の計算を行うが、各副領域の z 方向端部において参照する隣接格子点のデータは、もう一方の GPU 上のメモリ上に保存されている。そのため、副領域端部の計算において参照する隣接格子点のデータを他の GPU からあらかじめ袖領域(ghost domain)にコピーしておき、各副領域の計算終了後に袖領域を更新する方法を採用する。このとき、副領域とは別の連続なメモリを袖領域に割り当てることで、袖領域の更新を効率的に行う。
袖領域の更新には、
CUDA4.0 から可能となった Peer-to-Peer(P2P)通信を用いて、GPU 間で直接データ交換を行う。
本研究で行う3次元シミュレーションにおいては、各計算ノードに2基のCPU(IntelXeon E5-2690, 8cores, 2.9GHz)と2基のNVIDIA社のTESLA K20Xを備えたGPUクラスターを用いた。TESLA K20Xには，2688個のCUDAコアと6GBのビデオメモリが搭載されている。また、各ノードはInfiniband(QDR, 40Gbps)で通信可能となっている．計算プロ
グラムはCUDA Fortranを用いたため、コンパイラーにはPGI Accelerator Compiler v.12.10を用いた。
2.3 計量形態学に基づく 3 次元組織形態評価
実験的なミクロ組織観察で得られる組織形態の特徴を表す量（組織特徴値）としては、
2 次元観察では、粒径やその頻度分布、体積分率などが挙げられる。しかしながら、2次元観察では観察面の奥行き方向に存在する組織の形態は評価することができないた
め、真の組織特徴値を知ることは不可能である。そのため、緒言で述べたように、実験－ 84 －的なミクロ組織の 3 次元観察や計量形態評価が行われているところである。足立ら (4)は、機械研磨によるシリアルセクショニン
グと顕微鏡観察により得られるセクショニング画像（組織の 2 次元画像）をレンダリング（3 次元再構築）し、組織の 3 次元デジタル画像を得ることに成功している。さらに、得られた３次元デジタル画像を用いて、位相幾何学や微分幾何学に基づいて組織の形態評価を行い、組織の連結度や凹凸度などを定量的に、すなわち数値として評価することを可能としている。
一方、本研究で用いる MPF シミュレー

Fig.4 H-K plot.

ションでは、レンダリングなどの画像処理を行うことなく、計算結果としてミクロ組織形態の３次元デジタル画像を直接得ることができる。さらに、MPF シミュレーションを用いることで、ミクロ組織形態の時間変化も連続的な 3 次元デジタル画像として取得することができる。そこで本研究では、
MPF シミュレーションで得られる 3 次元デジタル画像を基に、足立らと同様の方法でミクロ組織表面のガウス曲率 K と平均曲率 H を算出し、ミクロ組織の形態およびその時間変化の定量評価を行う。
ミクロ組織表面の湾曲度を表すガウス曲率Kは、最大曲率1/r1と最小曲率1/r2の平均値として次式で与えられる。

K 11 1   2  r1 r2 (28)一方、平均曲率Hは、ミクロ組織表面の凹凸度を表す量であり、最大曲率と最小曲率の積として次式で与える。

H 1 1r1 r2(29)式(29)と式(30)で与えられるKとHを、ミクロ組織表面の全ての点について計算し、Fig.4に示すようなH-Kプロッ

Fig.5 Calculation of K and H using triangle patches.

トと呼ばれる分散図を作成することで、ミクロ組織形態がどのような形態かを定量的に評価できる。
H-K プロットにおいて、プロット点が第 1 象限の H2-K=0 の曲線と K=0 の間にある領－ 85 －

域 A にある場合には、組織表面が球形状(spheroid shape)または棒形状(cylindrical shape)の凸面であることを表している。また、原点近傍にプロットされる場合には、組織形態が直方体形状や板形状であることを示す。さらに、K<0 となる組織表面では、その組織形態が鞍形状(saddle shape)であることを示す。一方、H<0 となる場合にFig.6 Initial distribution of (a) austenite grains and (b)は、上記形状の凹面を表す。

austenite grain boundary.

Fig.5 に、MPF シミュレーショ
ンで得られるミクロ組織の 3 次元デジタル画像を基に H-K プロットを作成する方法を示す。MPF シミュレーションにより得られるi=0.5 の等値面を三角形パッチの集合体として表現し、STL(Standard Triangulated Language)ファイルで保存する。そして Fig.5 において A と示される、三角形パッチを構成する点 A における K と H の値は、点 A を囲む6 個の三角形パッチの面積 Si、各三角形パッチの内角i、各三角形パッチの法線が隣接する三角形パッチの法線となす角i を用いて、次式で計算する。
63 2  i i 1
K 

6Si 1H

6i

e 32i 16S4i 1(30)

i(31)

i

本研究では、MPF シミュレーションにより得られる 3 次元デジタル画像の三角形パッチ集合体への変換、K と H の算出を行うた
めに、株式会社サイバネッ
トシステムズの協力のもと MicroAVS ver.16.0 に曲率計算機能を追加した。

Fig.7 Time evolution of ferrite phase and carbon concentration
simulated by MPF simulation. Initial temperature is 1040K.
Cooling rate is 2.0K/s. Final temperature is 1020K. Vf represents
volume fraction of ferrite phase.

－ 86 －3. シミュレーション結果とK

考察
Fe-0.15wt.%C 合金の連続冷却過程における γ→α 変態の 3次元シミュレーションを行う。
Fig.6 に、母相であるオーステナイト相の初期多結晶構造を示す。解析領域の大きさは12.8×12.8×12.8m3 で あ り 、H128×128×128 の等間隔の差分格子で分割する。解析領域内

Kに 15 個のオーステナイト粒を配置し、古典的核形成理論で計算される速度で合計 41個のフェライト相を形成させる。初期温度は 1040K、冷却
速 度 は 2.0K/s 、 終 了 温 度 は1020K とする。また、界面エネルギーや界面易動度の異方H性は考慮しないため、フェライト相は等方的に成長する。

K

数値シミュレーションにおいては、全方向に周期境界条件を適用し、時間発展方程式の離散化には、空間 2 次精度の有限差分法、時間積分は 1 次精度のオイラー法を用いる。
Fig.7 に、MPF シミュレーシ
ョンで得られたγ→α変態にHよるフェライト相分布および炭素濃度の時間変化を示す。
また、フェライト相の冷却がFig.8 Contour plots of Gaussian curvature and Mean

進み温度が低下するとともに、 curvature and H-K plots for different volume fractions of
フェライト相が成長し、オー the ferrite phase.
ステナイト相においては、炭素濃度が増加していることが再現されている。
Fig.8 に、体積分率が 42%、64%、76%のフェライト相のガウス曲率 K および平均曲－ 87 －率 H を計算し、コンターマップとしてフェライト相表面に色付けした図とそれを基に作成した H-K プロットを示す。また、各 H-K プロットの下部と左側に、それぞれ H お
よび K の頻度分布を示す。フェライト相の体積分率が 42%の場合においては、H-K プ
ロットの第 1 象限に多くの点が存在する。これは、フェライト相はそれぞれ孤立し、球形状をしていることを表している。
次に、フェライト相の体積分率が 64%に増加すると、H-K プロットの第 3 象限のプロットが増加する。これはフェライト相が成長することで、孤立していたフェライト相が互いに連結し、負の曲率を持ったフェライト相の表面が形成されたためである。
さらに、フェライト相の体積分率が 76%まで増加すると、H-K プロットの第 3 象限のプロットがさらに増加するとともに、第 2 象限にもプロットが現れており、複数のフェライト相が接触し、連結度が増加していることが表されている。
以上の結果より、本研究では MPF シミュレーションで得られるミクロ組織形態の 3次元デジタル画像を用いて、ミクロ組織表面のガウス曲率 K や平均曲率 H を計算する方法を確立し、この方法で H-K プロットを作成することにより、ミクロ組織の連結度や凹凸度、さらにはその時間変化をも定量的に評価することができることが示された。
4. 結言
本研究では、Fe-C 合金で生じるγ→α変態挙動の高精度予測およびフェライト相形態変化の定量評価を実現するため、MPF 法による Fe-C 合金のγ→α変態シミュレーシ
ョンの並列 GPU 計算および MPF シミュレーションで得られるフェライト相の計量形態評価法を確立した。本研究により、MPF シミュレーションで得られるミクロ組織形態の 3 次元デジタル画像を用いて、ミクロ組織の連結度や凹凸度を表す、ガウス曲率 Kや平均曲率 H や H-K プロットの時系列変化を解析することが可能となった。MPF シミュレーションで得られるミクロ組織形態の 3 次元デジタル画像は、結晶塑性有限要素法などを用いた変形挙動解析の入力情報としても利用できるため、今後変形挙動解析と連携することにより、ミクロ組織形態の形成予測、3 次元組織形態の定評評価、組織形態と力学特性との関係を一貫して解析することが可能になるものと期待される。
謝辞
本研究は、2012 年度公益財団法人 JFE21 世紀財団技術研究助成鉄鋼技術研究により実施されたものであります。ここに深甚なる謝意を表します。株式会社サイバネットシステムズの松本陽司氏ならびに幸谷有子氏には、計量形態評価のために簡易可視化
ソフトウェア MicroAVS ver16.0 に曲率計算機能を追加して頂きました。厚く御礼申し上げます。

－ 88 －参考文献
(1) M. V. Kral, M. A. Mangan, G. Spanos and R. O. Rosenberg, Materials Characterization, 45,(2000), pp.17-23.
(2) R. Mendoza, I. Savin, K. Thornton and P. W. Voorhees, Nature Materials, 3, (2004),
pp.385-388.
(3) G. Spanos, D. J. Rowenhorst, A. C. Lewis and A. B. Geltmacher, MRS Bulletin, 23, (2008),
pp.597-602.
(4) 足立吉隆, ふぇらむ, 15 (2010), pp.28-35.
(5) Y-T. Wang, Y. Adachi, K. Nakajima and Y. Sugimoto, ISIJ International, 52, (2012),
pp.697-703.
(6) 足立吉隆, 軽金属, 61, (2011), pp.78-84.
(7) I. Steinbach and F. Pezzolla, A, Physica D, 134, (1999), pp.385.
(8) 高木知弘, 山中晃徳, フェーズフィールド法 –数値シミュレーションによる材料組
織設計-, 養賢堂, (2012).
(9) 小山敏幸, 高木知弘, フェーズフィールド法入門, 丸善出版, (2013).

－ 89 －