Meet up Chubu 「インフォマティクス」
ﾅｺﾞﾔｲﾉﾍﾞｰﾀｰｽﾞｶﾞﾚｰｼﾞ(ﾅﾃﾞｨｱﾊﾟｰｸ4F), 6/8(⽊), 2023.

フェーズフィールド法とデータ科学を活用した次世代材料設計小山敏幸
東海国立大学機構 名古屋大学 大学院工学研究科

1本日の内容 材料の設計とは？
 計算熱力学（ＣＡＬＰＨＡＤ法）
 計算組織学（フェーズフィールド法）
 イメージベース部材特性計算
 データ科学を活用した材料設計の加速2

材料設計の基盤的考え方
・第⼀原理電⼦状態計算
・分⼦動⼒学法(MD)・モンテカルロ法(MC)・CALPHAD法
・フェーズフィールド法・セルラーオートマトン
・有限要素法(FEM)・構成式

 空間および空間スケール
における複数階層
 マルチスケール
マルチフィジックス

原子構造(fcc, L12)

(γ+γʼ)組織

・数値流体⼒学(CFD), ・有限要素法(FEM)

タｰビンブレｰド[1] スケールブリッジング
 機械学習ジェットエンジン[2] データ同化[1] https://emira-t.jp/kengaku/5361/, (Oct., 2020).
[2] https://ja.wikipedia.org/wiki/ジェットエンジン, (Oct., 2020).

順問題・逆問題併⽤型の材料設計へ3【重要ポイント！】ｼﾞｪｯﾄｴﾝｼﾞﾝ内のﾀｰﾋﾞﾝ翼の内部組織
相:原子ｽｹｰﾙで均一な領域
組織:相の集合体
全ての素材の特性は、相と組織で決まる。
相の計算 ：CALPHAD法
組織の計算：ﾌｪｰｽﾞﾌｨｰﾙﾄﾞ法fcc(Ni)
L12(Ni3Al)

Ni-8.2at%Al-5.0at%Ti aged at 1023K
for (a)86.4ks and (b)345.6ks.
4材料設計とは？ ★材料特性の最適化（目的に合った材料の創製・高機能化）
方法論

・構成相の選択（発見も含む）
・構成相の分率の最適化
・組織サイズ（分布も含む）の最適化
・組織形態の最適化
5相と組織を制するもの、全ての材料開発のキャスティングボードを握る。

相の計算 ：ＣＡＬＰＨＡＤ法
組織の計算：フェーズフィールド法

いまや研究・開発で、世界的なディファクトスタンダードでありながら、大学で系統的に教えられていない。
（最近の20年で大きく進展した分野）
6

本日の内容 材料の設計とは？
 計算熱力学（ＣＡＬＰＨＡＤ法）
 計算組織学（フェーズフィールド法）
 イメージベース部材特性計算
 データ科学を活用した材料設計の加速

7

CALPHAD法：状態図計算（計算熱力学）
定番の状態図計算ソフト：(1) ThermoCalc(2) Pandat(3) CatCalc
Pandatのダウンロードサイト：
http://www.computherm.com/

入力ﾃﾞｰﾀ：構成相のギブスエネルギー関数（TDBファイル：ｷﾞﾌﾞｽｴﾈﾙｷﾞｰ関数の設定ﾌｧｲﾙ）→ 状態図計算、ｷﾞﾌﾞｽｴﾈﾙｷﾞｰ計算、・・・
TDBﾌｧｲﾙのｻﾝﾌﾟﾙ
NIMS MatNavi : 計算状態図ﾃﾞｰﾀﾍﾞｰｽ(CPDDB)
http://cpddb.nims.go.jp/8平衡状態図および準安定状態図を計算してみよう。

9

本日の内容 材料の設計とは？
 計算熱力学（ＣＡＬＰＨＡＤ法）
 計算組織学（フェーズフィールド法）
 イメージベース部材特性計算
 データ科学を活用した材料設計の加速10ﾌｪｰｽﾞﾌｨｰﾙﾄﾞ法 (現象論&連続体ﾓﾃﾞﾙ)
秩序変数の定義

ci (r, t ) : 保存変数(濃度,・・・)(0  ci  1)

si (r, t ) : 非保存変数(規則度,・・・) (  1  si  1)
全自由エネルギー汎関数

Gc {ci (r, t ), si (r, t )}Gsys     Egrad {ci (r , t ), si (r, t )}  Estr {ci (r, t ), si (r, t )} drr  E {c (r, t ), s (r, t )}  E {c (r, t ), s (r, t )}
imagii ele i発展方程式 Gsys
 ci ( r , t )    M ij , c j (r , t )tj Gsys si ( r , t )   Lij
 s j (r , t )
tj
11

12

当該分野の教科書・専門書・演習書

次世代材料デザイン手法の開発を目的に、材料学を横断する計算組織学分野を創る研究
13

演習書「ﾌｪｰｽﾞﾌｨｰﾙﾄﾞ法入門」内で解説されている組織形成

Phase decomposition

Martensitic transformation

Order-disorder phasetransition

Grain growth

Solidification
Grain boundary precipitation

Ferroic domain dynamics
14材料の内部組織形成を計算してみよう。

15

本日の内容 材料の設計とは？
 計算熱力学（ＣＡＬＰＨＡＤ法）
 計算組織学（フェーズフィールド法）
 イメージベース部材特性計算
 データ科学を活用した材料設計の加速16マクロ的力学特性の組織形態依存性
・次世代高強度材料開発では、組織形態情報の考慮が不可欠→組織形態情報（ﾓﾌｫﾛｼﾞｰ）の徹底活用の成否が、今後の高強度材料開発を著しく左右する。
[ISIJ研究会「計算工学による組織と特性予測技術II」（H22-24）]

12001000phase 1800600400

phase 00.05

1400
12001000phase 18006004000.100.150.2000.25

True strainc)Modulated structure
Both phases are deformed elastically
Phase 0 is deformed plasticallyBoth phases are deformed plastically

1400
12001000phase 1800600400

phase 02000

Both phases are deformed elastically
Phase 0 is deformed plasticallyBoth phases are deformed plasticallyphase 016000.050.100.15

True strain

0.200.250.000.05

SphereBoth phases are deformed elastically
Phase 0 is deformed plasticallyBoth phases are deformed plasticallyphase 11400
120010008006004000.100.150.2000.25

True straind)Layered structure

1600

Both phases are deformed elastically
Phase 0 is deformed plasticallyBoth phases are deformed plastically

1400
1200
1000

phase 18006004000a)phase 01600

200

phase 02000.00

1600

True stress / MPa

0.00

True stress / MPa

True stress / MPa

1600

Motteled structure

2002000b)

True stress / MPa

1400

1600

0.000.050.100.050.100.15

0.200.25

True strain

0.20

フェライト-パーライト鋼の応力-ひずみ曲線の組織依存性

Both phases are deformed elastically
Phase 0 is deformed plasticallyBoth phases are deformed plastically

12001000
80060040000.25

Modulated structure
Both phases are deformed elastically
Phase 0 is deformed plasticallyBoth phases are deformed plasticallyphase 11200
1000
800600
4000Motteled structurephase 11400

True strainc)

1400

0.15

phase 01600

phase 00.00

0.050.10d)

0.20

0.25

Layered structure
Both phases are deformed elastically
Phase 0 is deformed plasticallyBoth phases are deformed plasticallyphase 114000.15

True strain

12001000800600400

phase 02002000.00b)

200

True stress / MPa

SphereBoth phases are deformed elastically
Phase 0 is deformed plasticallyBoth phases are deformed plastically

True stress / MPaa)

True stress / MPa

True stress / MPa

1600

0.00

0.05

0.100.15

True strain

0.20

0.25

00.00

0.05

0.100.15

0.200.25

True strain

フェライト-ベイナイト鋼の応力-ひずみ曲線の組織依存性

17

時効析出型展伸用ﾏｸﾞﾈｼｳﾑ合金の組織・力学特性解析

Mg基hcp合金のﾌﾟﾘｽﾞﾏﾃｨｯｸGPｿﾞｰﾝ組織形
成PFｼﾐｭﾚｰｼｮﾝの実現。計算から得られる各時間の組織上で、転位が移動する最小応力（臨界分解せん断応力）を算出。
⇒ 析出硬化および過時効に対応する時効硬化曲線の計算。析出硬化期間は、ｿﾞｰﾝの体積分率増加期間とほぼ対応。
18

Fe-Nd-B希土類焼結磁石の磁気特性の組織依存性
ﾌｪｰｽﾞﾌｨｰﾙﾄﾞ法にて粒界相を含む多結晶不均一組織を計算し、その組織内の各位置に単相の磁気特性を割り付け、ﾏｲｸﾛﾏｸﾞﾈﾃｨｸｽに基づき、組織内の磁気ﾓｰﾒﾝﾄ場のﾀﾞｲﾅﾐｸｽおよび全体の磁気ﾋｽﾃﾘﾘｽをｲﾒｰｼﾞﾍﾞｰｽにて計算⇒ 微視的な磁化反転機構や磁気ﾋｽﾃﾘｼｽの組織形態依存性の解明が可能
以下の(a)～(c)は粒界相の磁性を、徐々に主相に近づけていった場合の計算結果(a) 10％(b) 30％2.02.00-1.02.0

aabcdefg1.0

Ms / T(c) 50％

bcdefg1.00-1.0h0-8.0-4.004.08.0g-1.0hh-2.0-2.0cdef1.0ab-2.0-8.0-4.00

4.0

H / MAm-1

8.0-8.0-4.00

4.0

8.0

19

PZTｾﾗﾐｸｽのMPB近傍における分極ﾄﾞﾒｲﾝ組織形成と誘電特性

組成c=0.5

MPB上におけるﾄﾞﾒｲﾝ組織変化

[10]
[01]

分極モーメントの向き

298Kにおける最大外部電場を[10]方向に30MVm-1で印加したときの平均組成0.5の分極ﾄﾞﾒｲﾝ組織変化ｼﾐｭﾚｰｼｮﾝ(a)t’ = 0, (b)t’ = 50, (c)t’ = 80, (d)t’ =250, (e)t’ = 500, (f)t’ = 540, (g)t’ = 555, (h) t’ = 750

PZTのﾓﾙﾌｫﾄﾛﾋﾟｯｸ相境界(MPB)における、外部電場下での分極ﾄﾞﾒｲﾝﾀﾞｲﾅﾐｸｽのﾌｪｰｽﾞﾌｨｰﾙﾄﾞ解析 ⇒ 電歪特性の組織および状態図依存性の解明が可能
20

ZrO2-4molY2O3における相分離と応力集中の組織依存性
ZrO2 - 8mol%YO1.5 at 1773Kb)548sa)219sa)c)1.64ksb)d)5.48kse)32.9ksf)110ks

(110)

100nm

[110]

[001]

[001]c)

(110)

100nm

[110]

Dark-field images of 8Y-PSZ aged at 1773 Kfor a)1800s, b)7200s and c)36ks.

ﾌｪｰｽﾞﾌｨｰﾙﾄﾞ法にて部分安定化ｼﾞﾙｺﾆｱの相分離組織形態を計算し、外力を作用させた際の内部応力場の組織形態依存性をﾌｪｰｽﾞﾌｨｰﾙﾄﾞ微視的弾性論にて計算⇒ 応力集中の組織形態（ｸﾗｯｸ等の欠陥も含む）依存性の解明が可能

21

PVME-PS系における相分離組織形成(森永研との共同研究)

Temperature/K523

Phase Diagram Calculated by Cell Model473b)

PVME:Mw=642.11k

PVME-40mass%PS at 413K (Calculation)a)1ks'b)2ksc)5ks'■ :Experimental(by Watanabe)■

PS:Mw=18.7k
■423

■

■

10μm

PVME-40mass%PS at 413K (Experimental)■ PS:Mw=44k■■a)15sb)30sc)60s

■

373

■

■ PS:Mw=280k:Binodal line(calc.):Spinodal line(calc.)

323
0PVMEa)

50

10μm

100
PSmass%PS
PVME-20mass%PS at 413K (Calculation)a)4ks'b)10ks'c)c)50ks

PVME-60mass%PS at 413K (Calculation)a)2.5ks'b)5ks'c)50ks

10μm

PVME-20mass%PS at 413K (Experimental)a)60sb)300s

10μm

PVME-60mass%PS at 413K (Experimental)c)600sa)300s

10μmb)900s

c)1.8ks

10μm22

本日の内容 材料の設計とは？
 計算熱力学（ＣＡＬＰＨＡＤ法）
 計算組織学（フェーズフィールド法）
 イメージベース部材特性計算
 データ科学を活用した材料設計の加速

23推定対象のパラメータ

不明な材料ﾊﾟﾗﾒｰﾀを、組織情報から推定!!

0仮想A-B-C三元系におけるスピノーダル分解を考える．

全て温度のみの関数とし，等温時効のシミュレーションデータを用いて推定する．

LA,C , 1LA,C , 0 LB,C , 1LB,C ,  , M 22

フィールド変数
濃度 ci (r, t )1：A成分，2：B成分，3：C成分．独立変数はc2，c3を定義（c1 = 1 – c2 – c3）

全自由エネルギー
Gsys  Gchem  Egrad ,（Gchem :化学的自由エネルギー， Egrad :勾配エネルギー）0101GccLccLL(cc)ccLL(cc)
RTclnc

ii
chem
12A,B
1
3A,CA,C
1
3
2
3B,CB,C23
 drr
i



ni
Li , j ：相互作用パラメータ222E 




(c)
dr
(c)
(c)
(c)(c)
dr


i
grad
2323r i 2r i ：勾配エネルギー係数

発展方程式M 22 Gsys Gsys ci    M ii 
 M ij  M ii 2  ichem    2 2 ci   2 c j    M ij  2   chem    2 ci  2 2 c j  j
 ci
 c j t

ただし， (i, j) = (2, 3), (3, 2)

M ij ：拡散易動度

 2chem  (1  2c2  c3 ) LA,B  c3 0 LA,C  c3 (2c2  3c3  2) 1LA,C  c3 0 LB,C  c3 (2c2  c3 ) 1LB,C  RT ln c2  ln(1  c2  c3 )
 chem022101
 3  c2 LA,B  (1  c2  2c3 ) LA,C  (6c3  6c2 c3  6c3  c2  2c2  1) LA,C  c2 LB,C  c2 (c2  2c3 ) LB,C  RT ln c3  ln(1  c2  c3 )

24

降下法[2]

計算方法［アジョイント法］【探索方向】準ニュートン法（BFGS公式）【ステップ幅】直線探索法（Armijo条件）

仮想A-B-C三元系におけるスピノーダル分解を考える．

状態ベクトル・シミュレーションモデル

評価関数

X(t )  (c2 (t ), c3 (t ),0 LA,C ,1LA,C ,0 LB,C ,1LB,C ,  , M 22 )T ,J X M (X)tM ; モデル演算子

tmax0 12

収束条件： Xk 1  Xk  1015

dt ,N meas.

3 N gridm 1i  2 n 1(n)(n)δ()()(t ) t

tct
c


_meas.
mii
22アジョイント方程式
λ
 M  λ   Xt

λ  tmax   0
λ  0 J
X  0 λ  t  ； アジョイント変数1) 濃度場
ci 2  chem
 2 ichem
j4  M ii   ci
 2 ci   M ij   ci  4 ci tci
ci
 2  chem
 2ichemj4
 M ji   c j 2 c j   M jj   c j  4 c j  
 ci
ci
ci



ただし， (i, j) = (2, 3), (3, 2)2) 相互作用パラメータ 0 LA,C ,1LA,C ,0 LB,C ,1LB,C
 chemLj2dr    M ij  cirLt
ij
ただし， L  0 LA,C ,1LA,C ,0 LB,C ,1LB,C3) 勾配エネルギー係数  2 M 22 c2 4 c2  M 22 c2 4 c3  M 23 c2 4 c2  2 M 23 c2 4 c3  dr
4444r
t2McMcM
c2M
c322323332333cccc33334) 拡散易動度 M 22
M 22t   c2 2  2chem   2c2  4 c2  c2  4 c3  drr[5] 矢部博，工学基礎 最適化とその応用，数理工学社（2016）

Y. Matsuura, Y. Tsukada, T. Koyama, Physical Review Materials, 5, 113801 (2021).

25計算条件［測定データ］
推定パラメータの真値

その他計算条件

LA,C  22000 (J/mol)
1LA,C  3000 (J/mol)
0LB,C  25000 (J/mol)
1LB,C 0 (J/mol)  5 1015 (Jm 2/mol)  1.0
M 22

温度
計算領域一辺の長さ計算領域の差分分割数
B, C成分の平均組成
相互作用パラメータ
拡散易動度
時間刻み t 

0900 K
100 nm1000.333, 0.333
LA,B  25000 (J/mol)  1.0, M ij  0.5
M 33
0.005

@900 K

PFシミュレーション結果

既知とする

測定データ

データ間隔 Δt’meas.

10データ数 nmeas.

5

測定時間 t’max

5026

推定結果
評価関数

相互作用パラメータ（0次項）相互作用パラメータ（1次項）

勾配エネルギー係数

拡散易動度

反復回数15程度で真値近傍に収束（反復回数30で真値に収束）

初期値

真値

LA,C  20000 → 22000 (J/mol)
LA,C 0 → -3000 (J/mol)0
LB,C  20000 → 25000 (J/mol)
1
LB,C 0 →0 (J/mol)
15
15
  6  10 →5  10 (J  m 2 /mol)  0.6 → 1.0
M 2201

また，上の初期値を他の値に変更した場合も推定が成功することを確認した．
27

熱拡散ｼﾐｭﾚｰｼｮﾝｼｽﾃﾑの開発

フェーズフィールド法による組織形態ﾋﾞｯｸﾞﾃﾞｰﾀの量産28

特徴量データの収集（組織形態およびS-S曲線）

ﾗﾒﾗの球状化２４５ 組織

7個の特徴量1. 第2相の体積分率
2. 界面積3. 種数 (穴の数)4. 連結度 (枝の数)5. 平均曲率Hの分散
6. ガウス曲率Kの分散
7. H2-Kの分散

２４５ S-S 曲線
1800

True stress / MPa

1600
1400

Both phases are deformed elastically
Phase 1 is deformed plasticallyBoth phases are deformed plasticallyMartensite
vf=0.91200

vf=0.71000

vf=0.5800

vf=0.3600

vf=0.1400

2個の特徴量1. 真応力2. 真ひずみ

Ferrite

200
00.00

0.05

0.10
0.15
True strain

0.20

0.25

これらのデータを、ニューラルネットワークに学習させる．
赤: 8個の入力データ , 青： 1個の出力データ

29

１４５個の組織情報と応力-ひずみ曲線のセットを用いてニューラルネットワーク学習
組織から７つの特徴量を抽出
Input layer

hidden layer

体積分率
界面積

Output layer

種数
連結度
平均曲率の分散ｶﾞｳｽ曲率の分散(H2-K)の分散

真ひずみ

真応力
出力は応力n 1iu m f   wij u nj   i  j 1uin : 第n層のニューロンiの出力値
wij : ニューロンjからニューロンiへの結合荷重

i : ニューロンiの閾値

入力層のニューロンの数 : 8
中間層のニューロンの数: 10
出力層のニューロンの数: 1

1
3層を持つ
f ( x) 1  exp(  ax )フィードフォーワード
ネットワーク30

ニューラルネット学習および感度解析結果
全入力因子の平均値を基準に、個々の因子を最小から最大まで変化させた時のSS曲線変化
体積分率

53.9%

連結度

界面積

種数

11.3%

平均曲率の分散

15.5%(H2-K)の分散

3.0%

2.4%

11.5%

ｶﾞｳｽ曲率の分散

2.4%

 体積分率、連結数、種数、界面積の順に応力-ひずみ曲線への寄与が大きい。
 体積分率は、高ひずみにおける応力に寄与、界面積、種数、連結度は、低ひずみでの応力分配に寄与。
 高強度鋼(ﾊｲﾃﾝ)における部材設計では、体積分率のみならず、ﾄﾎﾟﾛｼﾞｰﾊﾟﾗﾒｰﾀも重要な指標となる。
31

計算結果（フェーズフィールド法）

654個の組織においてそれぞれ3方向の熱伝導率を計算
計1962通り

緑： フィラー（断面は赤）32

各相の熱伝導率

計算結果（熱伝導率の計算）

フィラー： 0.02 Wm-1K -1

熱伝導率の異方性と組織形態の関係

図9 特徴的な部分の組織

母相： 0.12 Wm-1K -133

横長の結晶粒組織の場合（２次元計算）

34

DCGAN（深層畳み込み＋敵対的生成ﾈｯﾄﾜｰｸ）に基づく組織生成損失関数

本物を本物として正しく判断するための項

min max L( D ,  G )  Ex  pdata ( x )  ln D ( x; D )   Ez  pz ( z )  ln{1  D[G ( z; G ); D ]}GD

生成器から出力された画像を偽物として判断するための項

real images

判別器

D生成器

G
randomnoise

realorfake

fake images

TensorFlow 深い畳み込み生成的敵対的ネットワーク
https://www.tensorflow.org/tutorials/generative/dcgan

35

epoch=0epoch=50

epoch=100

epoch=150

epoch=200

epoch=300

DCGANによる生成器の学習過程

横長の結晶粒組織の場合（２次元計算）

36

DCGANで生成された画像を初期組織とし、PF法で短時間、追加計算

DCGANで生成された画像個々の結晶粒を抽出し、RGBで異なる色を塗る（255, 254, 253,・・・）。

PF法で短時間計算し、組織を均す。

GANで類似の組織が生成されるが、細部で物理的意味が消失。
⇒ PF法で短時間計算すると、界面張力のつり合いなど、リアルな組織へ改善!!

37次世代材料設計計算工学(a) 単相の物質ﾊﾟﾗﾒｰﾀ
・格子定数、弾性定数、界面ｴﾈﾙｷﾞｰ密度、ｷﾞﾌﾞｽｴﾈﾙｷﾞｰ、飽和磁化、ｷｭﾘｰ温度、拡散係数、・・・

第一原理計算, 原子・分子ｼﾐｭﾚｰｼｮﾝ,Calphad法, 各種ﾃﾞｰﾀﾍﾞｰｽ, ・・・[実験ﾃﾞｰﾀのﾃﾞｰﾀ同化]

組織ｼﾐｭﾚｰｼｮﾝ内の物質ﾊﾟﾗﾒｰﾀ

(b) 材料組織形態情報
・多結晶粒組織、析出組織、欠陥組織、転位組織、各種ﾄﾞﾒｲﾝ組織、・・・

組織情報

ﾌｪｰｽﾞﾌｨｰﾙﾄﾞ法, ｾﾙｵｰﾄﾏﾀｰ, LLG解析,
TDGL解析, ・・・[実験ﾃﾞｰﾀのﾃﾞｰﾀ同化]

組織の定量化[ ｻｲｽﾞ, 分率, 形状, 配向, 配置, ・・・(分布も含む)](c) 材料特性推定

構成式, FEM, FVM, BEM, 改良型ｾｶﾝﾄ法,均質化法, ﾏｲｸﾛﾏｸﾞﾈﾃｨｸｽ解析, ・・・・力学特性、磁気特性、電気特性、・・・[実験ﾃﾞｰﾀのﾃﾞｰﾀ同化]
ﾌﾟﾛｾｽｼﾐｭﾚｰｼｮﾝ内の材料ﾊﾟﾗﾒｰﾀ

(d) 部材・ﾃﾞﾊﾞｲｽの設計

38

まとめ
★計算熱力学（CALPHAD法）★計算組織学（フェーズフィールド法）★データ料学（機械学習・逆問題）

新たな体系化 通常のPCにて、⽇常的に計算できる環境
 多種・多量のデータの、効果的・相補的活⽤
 いつでも どこでも、順問題と逆問題
謝辞︓本研究の⼀部は、⽂部科学省ﾃﾞｰﾀ創出・活⽤型ﾏﾃﾘｱﾙ研究開発ﾌﾟﾛｼﾞｪｸﾄ事業（ﾃﾞｰﾀ創出・活⽤型磁性材料研究拠点）JPMXP1122715503 のご⽀援の下に⾏われたものであることを記し謝意を表します。
39