結晶粒成長過程のアダプティブ Phase-field シミュレーション
Adaptive Phase-Field Simulation during Grain Growth○正高木 知弘(神戸大・海科)

正冨田 佳宏(神戸大・工)Tomohiro TAKAKI, Faculty of Maritime Sciences, Kobe University, Kobe
Yoshihiro TOMITA, Graduate School of Science and Technology, Kobe University, Kobe
Eﬃcient adaptive ﬁnite element technique, in which ﬁne meshes are restricted around the grain boundaries, isapplied to phase-ﬁeld simulation of grain growth. The computing speeds for each step are improved with progressof simulation, because the number of ﬁnite elements is reduced with increasing of time steps due to the reductionof area for grain boundaries which is induced by the grain rotation and the movement of grain boundary.

Key Words : Grain Growth, Adaptive Finite Element Method, Phase-Field Method1

焼きなまし中の再結晶過程における微視組織の形態変化を把握することは，材料特性を高い精度で制御するために極めて重要である．結晶粒界のエネルギーを駆動力とした結晶粒の成長シミュレーションを可能とする
Phase-ﬁeld モデルが小林ら(1) によって提案されている．
Phase-ﬁeld 法は，不連続な界面を秩序変数の連続分布によって表し，界面の移動を秩序変数の時間発展として記述するため，複雑な形態変化を比較的容易に取り扱うことができる．しかしながら，界面における秩序変数の急峻な分布を許容するため，コンピュータシミュレーションにおいてかなり小さな計算格子を用いる必要があり，多くの計算コストを必要とする．本研究では，結晶粒成長過程を再現することのできる Phase-ﬁeld モデルにアダプティブ有限要素法(2) を適用し，数値シミュレーションの効率化を図り，その効果を検討する．

2

Phase-field モデル

本研究では，小林らによって提案された結晶粒成長過程を再現することのできる Phase-ﬁeld モデルを採用する(1) ．このモデルでは η と θ の 2 つの秩序変数を導入している．η は結晶粒内で 1，結晶粒界で 0 ≤ η < 1 の値を取る秩序変数．θ は各結晶粒の結晶方位であり結晶粒内で一定の値を取る．これらの秩序変数を用いた系の自由エネルギー汎関数は次のようになる．
 
F =f (η) +V

α2
ε2|∇η|2 + g (η) s |∇θ| +|∇θ|2 dV · · · · (1)22

ここで α，s，ε は正の定数である．f (η) はバルク中の自由エネルギー密度であり，結晶粒内 (η = 1) におい
2てのみ安定である f (η) = 12 (1 − η) を用いる．g (η) はg (0) = 0 である単調増加関数であり，g (η) = η 2 を用いる．η と θ は非保存量でり，各秩序変数の時間発展は次のようになる．
 2 2
∂η= Mη α ∇ η + (1 − η) − 2sη |∇θ| · · · · · · · · (2)
∂t


2 ∂θ
22 ∇θ
η
· · · · · · · · · · · · · · (3)= Mθ ∇ · ε ∇θ + sη∂t|∇θ|

ここで，Mη と Mθ は易動度である．
式 (2) と式 (3) の弱形式表示を行い，有限要素近似を行うと次に示す有限要素方程式が得られる．
 1   ∂η   2 Kη∂t

緒言+ Kη {η} = {fη } · · · · · · · · · · · · · · · · (4)Kθ1  ∂θ ∂t

+ Kθ2 {θ} = 0 · · · · · · · · · · · · · · · · · · · · (5)

ここで，二次元問題の場合，各マトリクスとベクトルは以下の通りである．
 1
1
Kη =





Kη2 =MηA 

α2

ϕN ϕM dA · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · (6)A{fη } =∂ϕN ∂ϕM∂ϕN ∂ϕM+∂x ∂x∂y ∂y+ (1 + 2s |∇θ|) ϕN ϕM dA (7)ϕM dA · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · (8)
A



1 2
η ϕN ϕM dA · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · (9)MθAKθ1 =

Kθ2 =ε2 +A

η2 s|∇θ|∂ϕN ∂ϕM∂ϕN ∂ϕM+∂x ∂x∂y ∂ydA · · · · · (10)時間に関してはクランク・ニコルソンの中央差分法を用いることで離散化を行う．この場合，時間 t + dt における秩序変数 η t+dt と θt+dt は次の連立一次方程式を解くことで得られる．
1  1
1  1
1  2   t+dt 1  2  tdt

Kη +1 dt1

Kθ +2Kη

2Kθ

1

3=

η  t+dt 2
θ

=dt

Kη −

1 dt


1

Kθ −2Kη

2Kθ

1

+{fη } (11)η  t2θ

· · (12)

アダプティブ有限要素法(2)要素の細分化および粗大化は図 1 のように，1 つの親要素を 4 つの子要素へ，4 つの子要素を 1 つの親要素へ変化させることにより行う．隣合う 2 つの要素の要素レベルが異なると，図 1 の節点 10 のように要素辺上に中間節点 (hanging node) が生じる．4 節点アイソパラ
メトリック要素を用いているため，節点 10 においては，η10 = 12 (η2 + η5 )，θ10 = 12 (θ2 + θ5 ) を満たす必要がある．ここでは，これらの条件式を隣合う要素のサイズの小さい
4
1
16 refinement of 4element 25221coarsening ofelements
132, 3, 4, and 5: node No.

95654hanging
node8111022373: element No.

Fig. 1 Reﬁnement and coarsening operations of element

Fig. 2 Initial mesh pattern and initial orientations

Fig. 4 Variations of the number of ﬁnite elements

要素（図 1 右図の要素 2 と 4）に対する式 (11) と式 (12)に代入することにより，中間節点の自由度を消去する．

4

シミュレーション

図 2 は本シミュレーションで用いた初期メッシュ分割と結晶粒の初期方位 (括弧内の値は方位差) を示している．
解析領域は，dX = dY = 3.2 とし，最も大きなレベル 0の要素で 10 × 10 分割している．レベル 0 からレベル 5 までの 6 段階の要素レベルを使用しているため，最も小さいレベル 5 の要素サイズは，dx = dy = 0.01 である．本
シミュレーションは以下の条件で行っている．α = 0.02，
s = 0.2，ε = 0.02，Mη = 10.，Mθ = 0.02，dt = 0.5．

(a)(b)(c)

Fig. 5 Variations of orientations for each grain

また，要素の細分化は要素内の 4 つの節点の η の平均値
ηave が ηave ≤ 0.8，または ηave > 0.8 かつ |∆θ| ≥ 0.1 の場合に行い，それ以外の要素は粗大化する．
図 3 は各ステップにおける要素分割，η 分布，結晶方位 θ 分布を示している．隣合う結晶粒の方位差が小さいほど結晶粒界の幅が狭く，粒界内における η の最小値が大きいことが分かる．また，結晶粒界付近のみ細かい要素が使用されており，要素のアダプティブ作業が効率的に行われていることが確認できる．
図 4 は要素数の時間変化を示している．全ての領域を最も細かなレベル 5 の要素で分割した場合，約 10 万個の要素を必要とするが，アダプティブ法を採用すると初期段階でその約 15%程度の要素数で解析が可能である．ス
テップ数が増加すると，結晶粒の回転および結晶粒界の移動に伴い結晶粒界の面積が減少するため，要素数が減少していることが分かる．このため，各ステップの計算時間はシミュレーション後期において特に短くなる．
図 5 は各結晶粒の方位の時間変化を示している．隣合う結晶粒の方位差が小さい場合は回転が生じ，結晶粒の粗大化が生じていることが分かる．また，各結晶粒の方位変化は連続的であり，要素のアダプティブ作業で用いた条件が妥当であるといえる．
考察および結言は省略．

参考文献
(d)
η 0.0
θ 0.0

0.5
1.5

1.0
3.0

Fig. 3 Time evolution of adaptive meshes, grain boundaries, and orientations. (a) 100, (b) 800, (c) 1600,
and (d) 4000 steps

(1)
(2)

(3)

Kobayashi, R., Warren, J. A., and Carter, W. C.,
Physica D, 140 (2000), 141 – 150.
Takaki, T., Fukuoka, T., and Tomita, Y., Proceedings
CD-ROM of the Sixth World Congress on Computational Mechanics, 310, 2004.
山本裕介, 高木知弘, 冨田佳宏, 機論, 70-699, A(2004),
1551-1557.