フェーズフィールド凝固シミュレーションへのアダプティブ有限要素法の適用
Phase-Field Simulation of Solidification using Adaptive Finite Element Method
高木 知弘（神戸大・海科） 冨田 佳宏（神戸大・工）Tomohiro TAKAKI, Kobe University
Yoshihiro TOMITA, Kobe University
FAX: 078-431-6286 E-mail: takaki@maritime.kobe-u.ac.jp
Phase-ﬁeld method is successfully applied to the study of the solidiﬁcation. Since the method mustallow the continuous of the phase-ﬁeld inside the narrow solid-liquid interface, very ﬁne meshes areneeded. In this study, the adaptive ﬁnite element method is introduced in order to conduct thephase-ﬁeld simulation eﬃciently. The degree of freedom of the hanging node, which occurs on theside of an element, is eliminated by conducting the matrix operation.

1.

緒言
凝固過程におけるデンドライト形成など，形態発展問題のシミュレーション手法としてフェーズフィールド法が注目されている．フェーズフィールド法は界面位置を明確に決める必要が無く，比較的容易に界面発展を記述することができる．一方で，非常に狭い界面領域内でフェーズフィールドの急峻な分布を許容する必要があるため，非常に小さな計算格子を必要とする．本研究では，アダプティブ有限要素法をフェーズフィールド法に適用することにより，凝固シミュレーションの効率化を図ることを目的としている．

10 においては，φ10 = 0.5(φ2 + φ5 ) を満たす必要がある．
ここでは，この条件式を隣合う要素のサイズの小さい要素に対する式 (1) に代入することにより，中間節点の自由度を消去する．例えば，図 1 右の要素 2 の式 (1) を， t+dt 
a a a a 
 φ2t+dt 111213142.とし，φ10 = 0.5(φ2 + φ5 ) を代入すると，
 t+dt  (a11 + a14 ) /2 a12 a13 a14 /2   φ2t+dt 

アダプティブ有限要素法
Warren ら1) によって提案された二元合金の凝固問題に対するフェーズフィールドモデルを採用する (詳細は文献 [1,2] を参照)．
空間は有限要素法，時間はクランク・ニコルソンの中央差分法を用いることで離散化を行うと，時間 t + dt におけるフェーズフィールド φt+dt は次の連立一次方程式を解くことにより得ることができる2) ．
    
1
11  2   t+dt 1  2   t11dt

Kφ +Kφ2

φ=dt

Kφ −2Kφ+{fφ }φ(1)   ここで，マトリクス Kφ1 と Kφ2 は，それぞれ熱伝導方程式でいう熱伝導・熱容量マトリクスに対応する．
要素の細分化および粗大化は図 1 のように，1 つの親要素を 4 つの子要素へ，4 つの子要素を 1 つの親要素へ変化させることにより行う．ある要素とその要素に隣接する要素のサイズレベルが異なると，図 1 の節点 10 のように要素辺上に中間節点 (hanging node) が生じる．4
節点アイソパラメトリック要素を用いているため，節点4
1
16 refinement of 4element 2522

1
coarsening ofelements
132, 3, 4, and 5: node No.

95654

hangingnode8111022: element No.

Fig. 1 Element reﬁnement and coasening

373

a21a31a41b=

11b21b31b41a22a32
a42b12b22b32b42b13b23b33b43(a21 + a24 ) /2
(a31 + a34 ) /2
(a41 + a44 ) /2

 (b11 + b14 ) /2 b12=(b21 + b24 ) /2
(b31 + b34 ) /2
(b41 + b44 ) /2

a23
a33
a43

b22b32b42

b14b24b34b44

a22
a32
a42

b13b23b33
b43

φ
7t+dt

a24
a34
a44φ φ11t+dt 10 t   φ2   f 

a23a33a43φt7φt11φt10+(2)2f7f11f10

φ7t+dta24 /2a34 /2a44 /2φ φ11t+dt  t 5  φ2   f b14 /22tb24 /2b34 /2b44 /2

φ7t

 φ11t
φ5+(3)

f7f11f10を得ることができる．同様の作業を要素 4 に対しても行うことで，中間節点 10 の自由度を消去することができる．連立一次方程式の解法としては，EBE 双共役勾配法を採用している．
凝固問題では移動する界面を追従するために，要素のアダプティブ作業をダイナミックに行う必要がある．ここでは，四分木データ構造を採用することにより，効率的な要素の細分化および粗大化作業を行う3) ．

3.

シミュレーション結果
図 2 は文献 [1] と同様の等温過飽和凝固デンドライト生成過程をシミュレーションの対象として得られた結果を示す1) ．解析領域左下に配置した結晶核から凝固が開始し，右方向と上方向へデンドライトが成長している．ここで，界面異方性として 4 回対称性を用いているため 1/4 領域を解析対象としている．使用したア
ダプティブ要素のレベル数は 8(レベル 0 – 7) としている．解析領域は 128×128 µm，最小要素 (レベル 7) の 1
辺 dx7 = dy 7 = 0.05 µm，最大要素 (レベル 0) の 1 辺dx0 = dy 0 = 27 dx = 6.4µm であり，全領域を最小要素のみで分割した場合，全要素数は 6,553,600 となる．要素(a) 2000 step(b) 4000 step(a) 32000(c) 6000 step(b) 57600(d) 64000(c) 62080(d) 8000 step
0.37
0.47
Concentration [at. frac.](e) 66560
Cut off domain to conductsimulation efficiently

Fig. 2 Time evolution for isothermal solidification
simulation (T = 1577.4K const., ∆ = 0.85,
c0 = 0.4 at. frac.)

G = 0.134 [K/µm]
V = 0.0025 [m/s]c0 = 0.4 at. frac.
288.0µm 115.2µm
0.37
0.47
Concentration [at. frac.](f) 110080 step

Fig. 4 Time evolution for directional solidification
simulation
initial transient andcompetitive growth[×10 ] 5

steady-state
growthFig. 3 Variations of the number of elemetns andCPU time versus the number of elementsの細分化は，要素内のフェーズフィールドの平均値 φaveが (0.001 < φave < 0.9) の場合に行い，それ以外の領域では粗大化を行っている．図 2(a) 中には 2000 ステップ時の要素分割図も示している．界面領域にレベル 7 の要素が配置されていることが分かる．図 3 は図 2 のシミュ
レーション過程における要素数変化，および CPU 時間(Pentium 4, 3GHz, 1CPU) – 要素数関係を示している．
シミュレーションが進むとデンドライトが成長し，界面の面積が増加するため，要素数は 2 次曲線的に増加している．また，CPU 時間と要素数はほぼ線形関係を示しており，各ステップにおける要素アダプティブ作業が非常に効率的に行われていることを示している．
図 4 は一方向凝固シミュレーションの結果を示している (詳細は文献 [3] 参照)．初期平滑界面を左から 1.6 µmの位置に配置し，それよりも左側を固相，右側を液相としている．要素の細分化は，(0.001 < φave < 0.9) または (| ∇cave |≥ 0.01 and φ < 0.001) の場合に行い，それ以外の領域では粗大化を行っている．シミュレーションが進むと平滑界面が乱れセル形態が発生し，セルまたはデンドライトの競争的な成長を経て，周期構造を有する定常成長状態に達することが分かる．図 4(e) にはデンド
ライトチップ先端近傍の要素分割を示している．界面領域および濃度勾配の大きい領域にそれぞれレベル 7 と 6の要素が配置されていることがわかる．図 5 は要素数とチップ先端の速度変化を示している．平滑界面状態では

The number of elements

50.008(e)(f)

40.006(d)30.0042 V = 0.0025 m/s
100(c)(b)

0.002(a)50000
Step No.

100000

Tip velocity [m/s]planar interface0Fig. 5 Variations of the number of elemetns andtip velocity

要素数は非常に少ないが，界面が乱れると要素数とチップ速度が急激に増加していることが分かる．また一方向凝固シミュレーションをより効率的に行うため，64000ステップから 2650 ステップ毎に解析領域の左端を解析対象から外している．このため，要素数がステップ状に変化している．また，チップ先端速度は定常時に引張り速度 V となっている．
参考文献

1) J. A. Warren, Boettinger, W. J., Acta Metall.
Mater., 43-2, pp.689 – 703, 1995.
2) 高木知弘, 福岡俊道, 冨田佳宏: 日本機械学会論文集
A 編, 70, 128-135, 2004.
3) T. Takaki, T. Fukuoka, Y. Tomita Proceedings CDROM of the Sixth World Congress on Computational Mechanics, 310, 2004.