Java による実践的科学技術プログラミング()―拡散相分離の Phasefield シミュレーション―小. は じ

め

に

山

敏

幸

全自由エネルギーの計算式

A B 2 元系の a 相の原子 1 モル当たりの化学的自由エネ

ルギーは，正則溶体近似に基づき，
今回 は，先 回の プログ ラム を修正 して， 拡散 相分 離のGac (c, T )＝°GaA(T )(1－c )＋°GaB(T )c＋LaA, B c(1－c )Phase
field シミュレーションプログラムを 2 種類作成する．

＋RT{c ln c＋(1－c ) ln (1－c)}

1 つは 1 次元濃度プロファイルの時間発展を計算するプログにて与えられる(7)．c が B 成分のモル分率である．組織形成(1)

ラムで，もう 1 つは 2 次元濃度場の時間発展プログラムで過程の計算であるので， c は時間 t および組織内の位置 r ＝ある．いずれも合金系は相分離型の仮想的な A B2 元系を(x, y, z)の関数として c(r, t )と表記するのが正しいが，ここ想定し，不規則相における 2 相分離の計算を行う．読者はでは式の煩雑さを避けるために c と記す．°
Gax ( T )は純成分

是非とも，ご自身のパソコンで計算を試みられたい．そしてX (＝ A, B )の自由エネルギーで，温度 T の関数である． R

拡散現象が目の前にてリアルタイムで進行していく様子を楽はガス定数で，右辺最後の項が原子の配置のエントロピーにしんでいただきたい．以下，まず計算手法について説明し，起因するエネルギー項である．右辺第 3 項が混合の過剰エ

その後，ソースコードとともにプログラムの解説を行う．

ンタルピー項で，係数の LaA, B は相互作用パラメータと呼ばれる． LaA, B は一般的には温度および組成の関数であるが，

. 計 算手法

.. AB 2 元系における a 相の相分離の計算式

本計算では定数と仮定し LaA, B＝ 25[kJ/mol]とおいた．本計
算は同一結晶構造内での等温時効における拡散相分解を対象と す るの で ， 自由 エ ネル ギ ー の基 準 の 0 と し て改 め て ，

相分離する相を a 相とする．本計算は，不規則 A B 2 元°GaA(T )(1－ c )＋°GaB (T )c を採用しても計算結果は影響を受系における a 相の相分離が対象であるので，秩序変数は濃けない．(拡散相分離の計算では，拡散ポテンシャルの空間

度場のみであり，したがって，本計算の Phasefield 法(1)(4)

勾配が溶質移動の駆動力になる．拡散ポテンシャルを計算すは，従来の CahnHilliard の非線形拡散方程式に基づく相分る際に自由エネルギーの組成に関する 1 次項は定数となり離の計算手法に等しい(5)(6)．ここでは独立な秩序変数として(温度一定)，さらに拡散ポテンシャルの空間勾配を計算する

B 成分のモル分率を採用する． Phase field 法では，まずこ際に，この定数は消えてしまう)．したがって，ここでは，の独立な秩序変数を用いて全自由エネルギーを汎関数形式に

議論を簡単にするために，自由エネルギーの基準について，定式化し，この全自由エネルギーが最も効率よく減少するよ°
GaA(T )＝°GaB(T )＝0 と設定する．

うに，その秩序変数の発展方程式(この場合は 2 元系の非線

次に界面エネルギーについて説明する．実際の界面エネル

形拡散方程式)が定義される．この発展方程式を数値計算す

ギーには，界面位置における平均場の化学的自由エネルギーることによって，秩序変数の時間および空間発展(すなわちも含まれるので，界面エネルギーと濃度勾配エネルギーは正組織形成過程)が算出される(1)．以下まず，使用する全自由

確には異なるが，ここでは(界面エネルギー)≡(濃度勾配エエネルギーの計算式について説明する．

ネルギー)として話を進める．計算式は，Esurf＝1
1
kc(;cA)2＋ kc(;cB)22
2

独 物質・材料研究機構新構造材料センター，組織・特性計算グループ 主幹研究員(〒3050047 つくば市千現 121)
Practical Java Programming in Materials Science and Engineering () ―Phase-field Simulation of Diffusional Phase Decomposition―;
Toshiyuki Koyama (National Institute for Materials Science, Structural Metals Center, Tsukuba)
Keywords: Jave application, Windows software, phase decomposition, spinodal decomposition, Phase-field method, diffusion equation
2009年 5 月11日受理

入 門 講 座＝て説明する(9)．まずエネルギーは RT にて無次元化する．ま1
1kc{;(1－c )}2＋ kc(;c)22
2

た 2 次元計算領域の 1 辺の長さを L ，差分計算の分割数を＝kc(;c)2(2)

N とすると，差分による空間分割セルの一辺の長さ b1 は，にて与えられる(1)(6)．スピノーダル分解の計算では濃度勾配

b1＝L/N である．この b1 を用いて距離を無次元する(差分計

エネルギー係数 kc は通常，方向に依存しない定数と仮定さ算を容易にするため)．また時間は b 21 / D にて無次元化されれる．kc の値は，原子間相互作用パラメータ Q，界面エネルる( D は拡散係数)．本計算では等温時効を想定しているのギー密度 gs ，界面幅 d1 ，および相互作用距離(8) d2 を用いで，易動度 Mc( c0, T )は定数 Mc となる．また拡散係数と易て，近似的に見積もることが可能である(9)．本計算では，kc

動度の関係式(後述)より，Mc の次元は[m2/s]/[J/mol]であ＝ 5.0 × 10

－15 [J ･ m2 / mol ]と置いた． k

c の単位については，

ここでは化学的自由エネルギーの単位[J/mol]との整合性か
ら[ J ･ m2/ mol ]を採用しているが，モル体積で割って[ J / m ]としている場合が多い点を記しておく．なお本講座はプログ
ラミングに重点を置いているので，議論を簡単にするためにる．以上より 2 次元(x, y )における非線形拡散方程式を，次元と合わせて表現すると，

(

2

(

&c
Mc RT &2( msys/RT) &2( msys/RT)＝＋D&(x/b1)2&( y/b1)2
t
& 2
b 1/D

( )

以上から，全自由エネルギーは，化学的自由エネルギーと界面エネルギー(濃度勾配エネルギー)の和として，空間積分

)

(10)

となる．ここで濃度勾配エネルギー定数 kc は b 21RT にて無における汎関数形式にて，f[LaA, Bc(1－c)＋RT{c ln c＋(1－c ) ln (1－c )}次元化される．
自己拡散係数 D
X(X＝A, B)と Mc(c, T )にはr＋kc(;c)2]dr(3)Mc(c, T )＝と表現される．Gsys の c による変分は，変分原理に基づき，

dGsys &Gac＝
－2kc;2cdc&c(4)にて与えられ(1)(6)，これは拡散ポテンシャルと呼ばれる．


][J/mmol] ( 9 )となり，無次元化した方程式は，

弾性歪エネルギーは考慮しない．

Gsys＝) [] [&c&2msys &2msys
1
m 2/ s
＝M c
＋ 2 
＝
2
&t
&x
&y
s
J/mol

発展方程式(非線形拡散方程式)の計算式

(

B

A

(11)

の関係がある(1)(6)．D
X は固溶体中における X 原子の自己拡散係数であるので，正確には D
X は温度および濃度(媒質の固溶体の濃度)の関数であるが，ここでは最も単純に，D
A＝D
B＝ D と置き，かつ組成 c については合金組成 c0 にて近似計算に用いる非線形拡散方程式(1)は，&cdGsys＝;･ Mc(c, T );
&t
dc

D
D
(1－c )＋c c(1－c )(RTRT )する．したがって，式( 11 )を変形することにより，式( 10))(5)右辺の係数部分が，McRT/D＝c0(1－c0)と導かれる．実際に相分解シミュレーションを行う際に，拡散係数の値は，無次にて与えられる． Mc ( c, T )は原子の拡散の易動度で，組成元化された計算時間 t′≡t/(b 21/D)を実時間 t に変換する時におよび温度の関数であるが(1)，合金組成 cのみ必要になる．無次元化された非線形拡散方程式は最終的0 および温度 T のに

関数 Mc(c0, T ) と近似する．式( 4 )を代入すると，[{}]＝;･ M (c , T )[(&&Gc );c]－2M (c , T )k ; c(&c&2( msys/RT) &2( msys/RT)＝ c 0( 1－ c 0)＋&t ′&(x/b1)2&( y / b 1) 2

&c
&G ac
＝;･ Mc(c0, T );－2kc;2c
&t
&c
c

a
c
2

0

c

0

c

)

(13)

であるので，無次元化された計算時間 t ′にて議論する限り

4

においては，拡散係数の値自体は特に必要ではない．なお本(6)を得る．ここで，DÃ ≡Mc(c0, T )(&2G ac /&c 2)および KÃ ≡ M c ( c 0,T )kc と置くと，CahnHilliard の非線形拡散方程式&c＝;･(DÃ ;c)－2K
Ã ;4c
&t

(7)

に 一致す ること がわか る． DÃ は 相互拡 散係数 に他な らない．なお実際の相分離シミュレーションでは，式( 4 )の拡
散ポテンシャル msys≡dGsys/dc を位置の関数 msys(r, t )として数値計算し，式( 5 )を&c
＝Mc;2msys&t(8)と変形して，差分法を用いて数値計算した方が効率的である．つまり組織内の任意の点における拡散ポテンシャル

msys ( r, t )を式( 4 )から求め， msys ( r, t )について直接，差分計算を行う手法である．
最後に 2 次元計算を例に，拡散方程式の無次元化につい

ま て り あ 第48巻 第10号(2009)Materia Japan

図1

計算に用いた AB 2 元系状態図．



計算では核形成成長過程も計算するので，実際の計算では0.07 および csp＝0.21, cBi＝0.93 および csp＝0.79 である．さ濃度揺らぎ項も乱数を用いて導入している．

て，プログラムの説明に入ろう．



AB 2 元系状態図

. プログラムの説明と計算例

計算に用いた A B 2 元系の状態図を図に示す．ここでは a 相のみを考慮している(液相は考慮していない)．実線
がバイノーダル線，および点線がスピノーダル線である．相基本的なプログラムの構造は，先回説明した F _ curve に分離の計算は，高温の単相領域で溶体化した状態を初期状態等しい．まず 1 次元濃度プロファイルの時間発展のプログとし，それを 2 相領域に持ちきたして等温時効した時の相ラムを以下に示す．プログラム名は， PD1D _ 001.java であ分離過程を対象とする．時効温度は 1000 K とする．この温る．プログラムの実行方法等は，先回の説明を参照していた度における平衡組成とスピノーダル組成はそれぞれ， cBi ＝

だきたい．さて先回同様，説明をプログラム内に記す．

↓//[プログラム(PD1D_001.java)] //[インポート文] 
import java.awt.
;import java.awt.event.
;import java.io.
;(1.0－c0);
Mc＝c0
c_flu＝0.1;//拡散の易動度//濃度ゆらぎ振幅の最大値//---- 時間 0 における初期濃度プロファイル設定--------------------------------------------------//乱数によって生成する場合
prog.ini_comp_field();//prog.datin();//ファイルから読み込む場合

public class PD1D_001 extends Frame{//[グローバル変数] 
//組織 1 辺の分割数
static int ND＝512;//濃度の分割数
static int nd＝ND;//濃度の分割数－1
static int ndm＝ND－1;//Window 全体の幅
static int width;//Window 全体の高さ
static int height;//描画領域の幅
static int xwidth;//描画領域の高さ
static int yheight;//Window の枠の幅(左右および下)
static int insetx;//Window の枠の幅(上)static int insety;//pstatic double PI＝3.141592;//ガス定数
static double RR＝8.3145;//組織内の濃度データ配列
static double[]ch＝new double[ND];//自由エネルギー曲線画面のグラフィックスオブジェクト
static Graphics g;//計算時間(カウント数)
static double time1;//温度(K)
static double temp;//合金組成(モル分率)
static double c0;//[コンストラクタ] 
public PD1D_001(){//描画画面の横と縦の長さ(ピクセル単位)
xwidth＝800; yheight＝200;//描画画面の枠の長さ
insetx＝4; insety＝30;//描画 Window 全体の横の長さ
width＝xwidth＋insetx
2;//描画 Window 全体の縦の長さ
height＝yheight＋insetx＋insety;//描画 Window のセット
setSize(width, height);//描画 Window の描画部の色を白に設定
setBackground(Color.white);//描画 Window を見えるようにする
setVisible(true);
addWindowListener(new WindowAdapter(){
public void windowClosing(WindowEvent e){System.exit(0);}//Window を閉じる時の操作(Window の右上×の設定)
});
}
//

[メインプログラム] 
public static void main(String[]args)throws Exception{//例外処理は行わない
PD1D_001 prog＝new PD1D_001();//PD1D_001のインスタンス prog を生成

int i;
int ip, im;double delt;
double al;
double time1max;
double[]ck＝new double[ND];
double[]ch2＝new double[ND];
double mu_chem, mu_surf;
double c1, c2;
double L0;
double kappa_c;
double Mc;
double c_flu;
double cddtt;//整数//整数(i＋1, i－1)//時間刻み(無次元)//計算領域の 1 辺の長さ//計算時間(最大カウント数)//拡散ポテンシャル//組織内の濃度データ予備配列//各ポテンシャル//濃度(A1, B2)//原子間相互作用パラメータ//濃度勾配エネルギー係数//易動度//濃度場の揺らぎの大きさ//濃度の増分

double b1;
double c2ip, c2im;
double sumc, dc;//差分ブロックサイズ//差分において c2 を中心に，その左右の濃度//濃度場の総和，平均組成からのずれ//---- 濃度プロファイルの時間発展の計算------------------------------------------------------------while(time1＜＝time1max){//---- 濃度プロファイルの表示-----------------------------------------------------------------------------if((((int)(time1)500)＝＝0)){prog.repaint();}//カウント数が500の倍数おきに濃度プロファイル描画//---- 濃度プロファイルの保存-----------------------------------------------------------------------------if((((int)(time1)1000)＝＝0)){ prog.datsave();}//カウント数が1000の倍数おきに濃度プロファイル保存//if(time1＝＝3000.0){prog.datsave();}//カウント数が3000の時に濃度プロファイル保存//---- 拡散ポテンシャルの計算-----------------------------------------------------------------------------for(i＝0; i＜＝ndm; i＋＋){
ip＝i＋1; im＝i－1;//周期的境界条件
if(i＝＝ndm){ip＝0;} if(i＝＝0){im＝ndm;}//位置 i における c1 と c2c2＝ch[i]; c1＝1.0－c2;//差分においてc2の左右の濃度
c2ip＝ch[ip]; c2im＝ch[im];(c1－c2)＋Math.log(c2)－Math.log(c1); //化学ポテンシャル差
mu_chem＝L0
(c2ip＋c2im－2.0
//濃度勾配のポテンシャル
mu_surf＝－2.0
kappa_cc2);//拡散ポテンシャル
ck[i]＝mu_chem＋mu_surf;}//---- 濃度場の時間変化(非線形拡散方程式の差分解陽解法)--------------------------------for(i＝0; i＜＝ndm; i＋＋){
ip＝i＋1; im＝i－1; if(i＝＝ndm){ip＝0;}if(i＝＝0){im＝ndm;} //周期的境界条件(ck[ip]＋ck[im]－2.0//非線形拡散方程式
cddtt＝Mc
ck[i]);//ch2[i]＝ch[i]＋cddtt
//濃度場の時間発展
delt;
(2.0
ch2[i]＝ch[i]＋(cddtt＋c_flu
Math.random()－1.0))
delt;//濃度場の時間発展(濃度揺らぎを考慮)}//
[濃度場の収支の補正] 
//数値計算であるので濃度場の収支の補正を行う(実際には毎ステップ行う必要はない) 
sumc＝0.; for(i＝0; i＜＝ndm; i＋＋){sumc＋＝ch2[i];} //濃度プロファイルの積分//濃度プロファイルの変動量
dc＝sumc/(double)nd－c0;
for(i＝0; i＜＝ndm; i＋＋){ch[i]＝ch2[i]－dc;
if(ch[i]＞＝1.){ch[i]＝1.0－1.0e－6;}if(ch[i]＜＝0.){ch[i]＝1.0e－6;}}//濃度場の補正//濃度が1を超えた場合の補正//濃度が0を切った場合の補正//
[時間増加] 
time1＝time1＋1.0;}//while//計算時間の増加// ---------------------------------------------------------------------------------------------------System.out.printf(``¥n 終了しました．グラフの右上×をクリックして終了してください．¥n'');}//main//---- 各種パラメータ設定------------------------------------------------------------------------------------//[K]temp＝1000.0;//合金の平均組成(ここでは A
c0＝0.4;40 atB 合金を設定している)
delt＝0.01;//計算の繰返し回数
time1＝0.0;//繰返し回数の最大値
time1max＝1.0e＋08;//以下はサブルーチン//初期濃度場設定 
public void ini_comp_field(){
int i;
double rnd0, fac1;//初期濃度揺らぎの最大変化量を1に設定
fac1＝0.01;
for(i＝0; i＜＝ndm; i＋＋){
rnd0＝2.0
Math.random()－1.0; ch[i]＝c0＋rnd0
fac1; //乱数にて初期濃度場を設定}}

al＝500.0;al＝al
1.0e－9;b1＝al/(double)nd;//2 次元計算領域の 1 辺の長さ(nm)//(m)に変換//差分 1 ブロックのサイズ

L0＝2.5e＋04;
L0＝L0/RR/temp;//原子間相互作用パラメータ(J/mol)//無次元化//濃度勾配エネルギー係数，単位は[J m^2/mol]
kappa_c＝5.0e－15;
kappa_c＝kappa_c/b1/b1/RR/temp; //(b1^2
rrtemp)で無次元化






//

濃度場描画 
public void paint(Graphics g){
int i, ii, i1, ii1, i2, ii2;
double x, xmax, xmin, dx, y, ymax, ymin, dy, d0;
double c, cmax, cmin;double gx1, gy1, gx2, gy2;
int ixmin, iymin, ixmax, iymax, igx1, igy1, igx2, igy2, id0;
int idx, idy;
double col;



入 門 講 座//横軸の最小値，最大値，分割間隔(実際の値)
xmin＝0.; xmax＝1.; dx＝0.1;//縦軸の最小値，最大値，分割間隔(実際の値)
ymin＝0.; ymax＝1.; dy＝0.1;//濃度の最小値，最大値(実際の値)
cmin＝0.0; cmax＝1.0;
ixmin＝0; iymin＝0; ixmax＝xwidth;iymax＝yheight;//ピクセル単位の場合
idx＝(int)(0.1
ixmax);
idy＝(int)(0.1
iymax);

g.drawLine(insetx＋igx1, insety＋igy1, insetx＋igx2, insety＋igy2);//濃度プロファイルの描画
}
}
//

[データの保存] 
private void datsave()throws Exception{
int i;//色を白に設定
g.setColor(Color.white);g.fillRect(insetx, insety, xwidth, yheight); //画面を上で指定した色で塗る//色を黒に設定
g.setColor(Color.black);g.drawRect(insetx, insety, xwidth, yheight); //グラフの外周を上で指定した色で描く//グラフ内の縦横の線を等間隔に描く
for(i＝0; i＜＝ixmax; i＋＝idx){g.drawLine(insetx＋i, insety＋iymin,
insetx＋i, insety＋iymax);}for(i＝0; i＜＝iymax; i＋＝idy){g.drawLine(insetx＋ixmin, insety＋i,
insetx＋ixmax, insety＋i);}//差分ブロックの長さの半分
d0＝1.0/(double)nd/2.0;
id0＝1＋(int)(((double)ixmax－(double)ixmin)/(xmax－xmin)
d0);//d0のピクセル化
System.out.printf(``f ¥n'', time1);//保存ファイル名を test.dat とする．
PrintWriter outfile＝new PrintWriter(
new BufferedWriter(new FileWriter(``test.dat'', true)));}//[データの読込み] 
private void datin()throws Exception{
int i;
String s_data;//計算の繰返し回数を標準入出力に表示//色を赤に設定
g.setColor(Color.red);for(i＝0; i＜＝ndm; i＋＋){//濃度プロファイルの値(実際の値)
i1＝i; i2＝i＋1;(double)i1＋d0; gx2＝1./(double)nd
(double)i2＋d0;
gx1＝1./(double)nd
ii1＝i1; ii2＝i2; if(i＝＝ndm){ii2＝0;}
gy1＝(ch[ii1]－cmin)/(cmax－cmin); gy2＝(ch[ii2]－cmin)/(cmax－cmin);//濃度プロファイルの値をスクリーン座標に変換(gx1－xmin)/(xmax－xmin)
igx1＝(int)(((double)ixmax－(double)ixmin)
＋(double)ixmin);igy1＝(int)((double)iymin＋(double)iymax－(((double)iymax－(double)iymin)/(ymax－ymin)
(gy1－ymin)＋(double)iymin));(gx2－xmin)/(xmax－xmin)
igx2＝(int)(((double)ixmax－(double)ixmin)
＋(double)ixmin);igy2＝(int)((double)iymin＋(double)iymax－(((double)iymax－(double)iymin)/(ymax－ymin)
(gy2－ymin)＋(double)iymin));//ファイルのオープン(追記)//カウントの書き込み
outfile.println(time1);
for(i＝0; i＜＝ndm; i＋＋){outfile.println(ch[i]);}//濃度場の書き込み//ファイルのクローズ
outfile.close();

BufferedReader infile＝new BufferedReader(new FileReader(``ini000.dat''));//ファイルのオープン
s_data＝infile.readLine();time1＝new Double(s_data).doubleValue();for(i＝0; i＜＝ndm; i＋＋){s_data＝infile.readLine();ch[i]＝new Double(s_data).doubleValue();}
infile.close();//文字列として読み込み//文字を数値へ変換//文字列をして読み込み//文字を数値へ変換//ファイルのクローズ}// 
}//PD1D_001//プログラム終了 



本計算では，初期濃度プロファイルを乱数にて設定してい長型相分離の特徴を示すようになる．最後に合金組成がバイるので，最後のサブルーチンは使用していない．図はプロ

ノーダル線のすぐ内側の c0＝0.1では，完全に核形成成長型

グラム実行結果の 1 例で， A 40 at  B 合金の 1000 K 等温

相分離の相分離となる．通常，スピノーダル分解と核形成

時効における濃度プロファイルの時間発展である．初期に典成長型相分解は，異なる相分離メカニズムとして説明される型的なスピノーダル分解が生じ，時効の進行に伴い，濃度ピ場合が多いが，原子の拡散の観点からは，基本的には 1 つーク同士が互いに溶質を奪い合ってオストワルド成長していの非線形拡散方程式から計算される単一の現象である(5)．たく様子が計算されている．

だし，図 3 の中央組成の濃度プロファイルと端の組成の濃次に合金組成を変えて，同様に相分解させた時の濃度プロ度プロファイルの比較からわかるように，得られる組織形態

ファイル(時間 t ′＝ 3000 における結果)を図に示す．状態

は，スピノーダル線付近を境に大きく変化する．このこと

図の中央組成付近では，典型的なスピノーダル分解組織となが，従来，スピノーダル分解と核形成成長型相分解を分けるが，合金組成が低下して，ほぼスピノーダル線上の c0 ＝て考えてきた所以であろう(もちろんスピノーダル点は，過0.2 では，濃度プロファイルの周期性が若干崩れてくる．合

飽和固溶体からの微小濃度揺らぎに対して，相分離に対する金組成がスピノーダル線とバイノーダル線の中央付近である

化学的駆動力が負から正に変わる点であり，相分離の本質的

c0 ＝ 0.15 では，まばらに濃度ピークが形成され，核形成成

な機構は異なる)．

図2

A 40 at  B 合金の 1000 K 等温時効における 1
次元濃度プロファイルの時間変化．
＝10, (b) t′＝20, (c) t′＝50, (d) t′＝200,(a ) t ′(e ) t ′＝1000, (f) t′＝3000.

ま て り あ 第48巻 第10号(2009)Materia Japan

図3＝3000)におけるAB 合金の 1000 K 等温時効(t′1 次元濃度プロファイルの合金組成による変化．
(a) A50 atB, (b) A40 atB, (c) A30 atB,(d) A20 atB, (e) A15 atB, (f) A10 atB.

さて次に，このプログラムを 2 次元に拡張してみよう．
説 明 は 以 下 の よ う に な る ( プ ロ グ ラ ム 名 は PD2D _001.java )．先の濃度プロファイルの計算プログラムと比較しながら見ていただきたい．

↓//[プログラム (PD2D_001.java)] //[インポート文] 
import java.awt.
;import java.awt.event.
;import java.io.
;public class PD2D_001 extends Frame{//[グローバル変数] 
//組織 1 辺の分割数
static int ND＝64;//濃度の分割数
static int nd＝ND;//濃度の分割数－1
static int ndm＝ND－1;//Window 全体の幅
static int width;//Window 全体の高さ
static int height;//描画領域の幅
static int xwidth;//描画領域の高さ
static int yheight;//Window の枠の幅(左右および下)
static int insetx;//Window の枠の幅(上)static int insety;//pstatic double PI＝3.141592;//ガス定数
static double RR＝8.3145;//組織内の濃度データ配列
static double[][]ch＝new double[ND][ND];//自由エネルギー曲線画面のグラフィックスオブジェクト
static Graphics g;//計算時間(カウント数)
static double time1;//温度(K)
static double temp;//合金組成(モル分率)
static double c0;//[コンストラクタ] 
public PD2D_001(){//描画画面の横と縦の長さ(ピクセル単位)
xwidth＝400; yheight＝400;//描画画面の枠の長さ
insetx＝4; insety＝30;//描画 Window 全体の横の長さ
width＝xwidth＋insetx
2;//描画 Window 全体の縦の長さ
height＝yheight＋insetx＋insety;//描画 Window のセット
setSize(width, height);//描画 Window の描画部の色を白に設定
setBackground(Color.white);//描画 Window を見えるようにする
setVisible(true);
addWindowListener(new WindowAdapter(){public void windowClosing(WindowEvent e){System.exit(0);}//Window を閉じる時の操作(Windowの右上×の設定)
});
}
//

[main プログラム] 
public static void main(String[] args) throws Exception{//例外処理は行わない
PD2D_001 prog＝new PD2D_001();//PD2D_001のインスタンス prog を生成

int i, j;
int ip, im, jp, jm;
double delt;
double al;
double time1max;
double[][]ck＝new double[ND][ND];
double[][]ch2＝new double[ND][ND];
double mu_chem, mu_surf;
double c1, c2;
double L0;
double kappa_c;
double Mc;
double c_flu;
double cddtt;//整数//整数(i＋1, i－1, j＋1, j－1)//時間刻み(無次元)//計算領域の 1 辺の長さ//計算時間(最大カウント数)//拡散ポテンシャル//組織内の濃度データ予備配列//各ポテンシャル//濃度(Fe1, Cu2)//原子間相互作用パラメータ//濃度勾配エネルギー係数//易動度関数とその微分//濃度場の揺らぎの大きさ//濃度の増分//---- 濃度場の保存----------------------------------------------------------------------------------------------if((((int)(time1)500)＝＝0)){prog.datsave();}//カウント数が500の倍数おきに濃度場保存//if(time1＝＝3000.0){prog.datsave();}//カウント数が3000の時に濃度場保存//---- 拡散ポテンシャルの計算-----------------------------------------------------------------------------for(i＝0; i＜＝ndm; i＋＋){for(j＝0; j＜＝ndm; j＋＋){
ip＝i＋1; im＝i－1; jp＝j＋1; jm＝j－1;//周期的境界条件
if(i＝＝ndm){ip＝0;} if(i＝＝0){im＝ndm;}//周期的境界条件
if(j＝＝ndm){jp＝0;} if(j＝＝0){jm＝ndm;}//位置(i, j)における c1 と c2c2＝ch[i][j]; c1＝1.0－c2;c2ip＝ch[ip][j]; c2im＝ch[im][j];c2jp＝ch[i][jp]; c2jm＝ch[i][jm];//差分において c2 の前後左右の濃度(c1－c2)＋Math.log(c2)－Math.log(c1);mu_chem＝L0
//化学ポテンシャル差(c2ip＋c2im＋c2jp＋c2jm－4.0
mu_surf＝－2.0
kappa_cc2);//濃度勾配のポテンシャル//拡散ポテンシャル
ck[i][j]＝mu_chem＋mu_surf;}}//---- 濃度場の時間変化(非線形拡散方程式の差分解陽解法)--------------------------------for(i＝0; i＜＝ndm; i＋＋){for(j＝0; j＜＝ndm; j＋＋){
ip＝i＋1; im＝i－1; jp＝j＋1; jm＝j－1;//周期的境界条件
if(i＝＝ndm){ip＝0;} if(i＝＝0){im＝ndm;}//周期的境界条件
if(j＝＝ndm){jp＝0;} if(j＝＝0){jm＝ndm;}(ck[ip][j]＋ck[im][j]＋ck[i][jp]＋ck[i][jm]－4.0ck[i][j]);cddtt＝Mc
//非線形拡散方程式//ch2[i][j]＝ch[i][j]＋cddtt
//濃度場の時間発展
delt;
(2.0
ch2[i][j]＝ch[i][j]＋(cddtt＋c_flu
Math.random()－1.0)) 
delt;//濃度場の時間発展(濃度揺らぎ導入)
}
}
//

[濃度場の収支の補正] 
//数値計算であるので濃度場の収支の補正を行う(実際には毎ステップ行う必要はない)．
］ 
sumc＝0.;for(i＝0; i＜＝ndm; i＋＋){for(j＝0; j＜＝ndm; j＋＋){//濃度場の積分
sumc＝sumc＋ch2[i][j];}}//濃度場の変動量
dc＝sumc/(double)nd/(double)nd－c0;
for(i＝0; i＜＝ndm;i＋＋){for(j＝0; j＜＝ndm;j＋＋){ch[i][j]＝ch2[i][j]－dc;
if(ch[i][j]＞＝1.){ch[i][j]＝1.0－1.0e－6;}if(ch[i][j]＜＝0.){ch[i][j]＝1.0e－6;}}}//濃度場の補正//濃度が 1 を超えた場合の補正//濃度が 0 を切った場合の補正//
[時間増加] 
//計算時間の増加
time1＝time1＋1.0;}//while//---------------------------------------------------------------------------------------------------------------------------------------------------//差分ブロックサイズ
double b1;
double c2ip, c2im, c2jp, c2jm;//差分ブロックにおいて c2 を中心に，その上下左右の濃度//濃度場の総和，平均組成からのずれ
double sumc, dc;//---- 各種パラメータ設定-----------------------------------------------------------------------------------//[K]temp＝1000.0;//合金の平均組成(ここでは A
c0＝0.4;40 atB 合金を設定している)//時間きざみ
delt＝0.04;//計算の繰返し回数
time1＝0.0;//繰返し回数の最大値
time1max＝1.0e＋08;
al＝60.0;al＝al
1.0e－9;b1＝al/(double)nd;//2 次元計算領域の 1 辺の長さ(nm)//(m)に変換//差分 1 ブロックのサイズ

L0＝2.5e＋04;
L0＝L0/RR/temp;//原子間相互作用パラメータ(J/mol)//無次元化//以下はサブルーチンである．
//
初期濃度場設定 
public void ini_comp_field(){
int i, j;double rnd0, fac1;//初期濃度揺らぎの最大変化量を 1に設定//乱数にて初期濃度場を設定}//拡散の易動度// 
//チラツキ防止のためdraw関数を別に定義
public void update_draw(Graphics g){ g＝getGraphics(); paint(g);}//---- 時間 0 における初期濃度場設定------------------------------------------------------------------prog.ini_comp_field();//prog.datin();//ファイルから読み込む場合//
濃度場描画 
public void paint(Graphics g){//g.clearRect(0, 0, width, height);//Window をクリア//---- 濃度場の時間発展の計算-----------------------------------------------------------------------------while(time1＜＝time1max){//---- 濃度場の表示----------------------------------------------------------------------------------------------//カウント数が200の倍数おきに濃度場描画
if((((int)(time1)200)＝＝0)){prog.update_draw(g);}//描画時のチラツキを抑えるため//if((((int)(time1)200)＝＝0)){prog.repaint();}






}//main

fac1＝0.01;
for(i＝0; i＜＝ndm;i＋＋){for(j＝0; j＜＝ndm;j＋＋){
rnd0＝2.0
Math.random()－1.0;ch[i][j]＝c0＋rnd0
fac1;
}
}

//濃度勾配エネルギー係数，単位は[J m^2/mol]
kappa_c＝5.0e－15;
kappa_c＝kappa_c/b1/b1/RR/temp;//濃度勾配エネルギー係数，(b1^2
rrtemp)で無次元化(1.0
Mc＝c0
c0);
c_flu＝0.1;

System.out.printf(``¥n 終了しました．グラフの右上×をクリックして終了してください．¥n'');int i, j, ii, jj;
double c, x, xmax, xmin, y, ymax, ymin, rad0;
int ixmin＝0, iymin＝0, igx, igy, irad0;
int ixmax＝xwidth, iymax＝yheight;
int icol;xmin＝0.; xmax＝1.;
ymin＝0.; ymax＝1.;//横軸の最小値，最大値//縦軸の最小値，最大値//差分ブロックの長さの半分
rad0＝1.0/(double)nd/2.0;
irad0＝1＋(int)(((double)ixmax－(double)ixmin)/(xmax－xmin)
rad0);//rad0 のピクセル化

入 門 講 座System.out.printf(``f ¥n'', time1);//計算の繰返し回数を標準入出力に表示

for(i＝0; i＜＝nd; i＋＋){for(j＝0; j＜＝nd; j＋＋){//濃度場の位置座標(実際の値)(double)i＋rad0;
x＝1.0/(double)nd
(double)j＋rad0;
y＝1.0/(double)nd//濃度場の位置座標(スクリーン座標に変換)(x－xmin)/(xmax－xmin)
igx＝(int)(((double)ixmax－(double)ixmin)
＋(double)ixmin);(y－ymin)/(ymax－ymin)igy＝(int)(((double)iymax－(double)iymin)
＋(double)iymin);//カウントの書き込み

outfile.println(time1);
for(i＝0; i＜＝ndm;i＋＋){for(j＝0; j＜＝ndm;j＋＋){
outfile.println(ch[i][j]);}}outfile.close();//濃度場の書き込み//ファイルのクローズ}//[データの読込み] 
private void datin()throws Exception{
int i, j;String s_data;//個々の差分ブロックの濃度値
ii＝i; jj＝j;//周期的境界条件
if(i＝＝nd){ii＝0;} if(j＝＝nd){jj＝0;}(1.0－ch[ii][jj]));//色の諧調をグレースケールにする
icol＝(int)(255.0
//icol＝(int)(255.0
//明暗を反転する場合
ch[ii][jj]);
if(icol＞＝255){icol＝255;}if(icol＜＝0){icol＝0;}//明暗の範囲の補正//濃度を明暗で設定
g.setColor(new Color(icol, icol, icol));g.fillRect(insetx＋igx－irad0, insety＋igy－irad0, irad0
2, irad02);//個々の差分ブロック描画

BufferedReader infile＝new BufferedReader(new FileReader(``ini000.dat''));//ファイルのオープン
s_data＝infile.readLine();time1＝new Double(s_data).doubleValue();for(i＝0; i＜＝ndm; i＋＋){for(j＝0; j＜＝ndm; j＋＋){s_data＝infile.readLine();ch[i][j]＝new Double(s_data).doubleValue();}}
infile.close();}}}//文字列をして読み込み//文字を数値へ変換//文字列をして読み込み//文字を数値へ変換//ファイルのクローズ}//[データの保存] 
private void datsave() throws Exception{
int i, j;//保存ファイル名を test.dat とする．
PrintWriter outfile＝ new PrintWriter(new BufferedWriter(new FileWriter(``test.dat'', true)));//ファイルのオープン追記// 
}//PD2D_001//プログラム終了 



このプログラムを実行すると，2 次元濃度場の時間発展が読み出して，濃度プロファイルや濃度場を表示するプログラ

リアルタイムにてパソコン画面で見ることが出来る．図にム，個々の時間の濃度プロファイルデータを個別のファイル

計算結果の濃度場の時間発展を示す．図 4(a)～(h)は，合金に取り出すプログラム(図 2 や図 3 のような図を書く時に必

組成 c0 ＝ 0.4 ( A 40 at  B )の時効温度 1000 K におけるスピ要)，また濃度場の 2 次元数値データを jpg 形式などのイメ

ノーダル分解の時間発展である．図中の相分解組織の明暗がージファイルに変換して保存するプログラムなども必要であ

局所的な B 成分の組成に対応し，白～黒が 0 at  B ～ 100

る．これらのプログラムに関しては，紙面の関係で本講座でat  B を表している．また計算領域 1 辺は 60 nm に設定しは説明しないが，これら作業用のプログラムも，本掲載プロてある．

グ ラ ム の ダ ウ ン ロ ー ド ペ ー ジ ( http: // www.nims.go.jp /初期状態は過飽和固溶体(最大で± 1 程度の濃度揺らぎ

mpsg/PhaseField_Modeling.htm )に公開するので，活用しを乱数によって与えてある)であり，相分解の初期に均一なていただけると嬉しい(ただし公開しているプログラム等に“まだら構造”が形成され[( a )( d )]，その後，組織は孤立

関して不具合・トラブルがあっても，著者および社団法人日した微細な析出粒子が分散する組織となり[(e)(g)]，時効

本金属学会は責任を負いませんのでご了承ください)．

の進行に伴いオストワルド成長によって析出相が粗大化していく[(g)(h)]．合金組成や，時効温度を変えて種々の条件

. お

わ

り

に

下での相分解過程を，パソコン上にてその場観察しながら計算を進めることができる．

以上，AB2 元系の拡散相分離シミュレーションを例にとさて以上は，相分離過程を計算するプログラムであるが，り，パソコン上で簡潔に相分解シミュレーションを行うプロ

実際に各種作業を進めるためには，保存された計算データを

グラムについて説明した．本プログラムは，学部や大学院における材料組織学の授業においても効果的に活用できると思う．最近では，多くの授業で PowerPoint による説明が行われるので，その場で，濃度プロファイルや濃度場の時間発展をデモンストレーションすることは教育効果としても有効であろう．もちろん，本プログラムを少し改良するだけで，現
実 の 2 成 分 系相 分 離 のプ ロ グ ラム を 作 成す る こと が で きる．例えば， Fe Cr や Fe Cu のフェライト相のスピノーダ
ル分解プログラム(9)は比較的容易に作成可能である(ただしGibbs エ ネ ル ギ ー 内 の 磁 気 項 の 取 り 扱 い は 若 干 複 雑 で ある)．また本計算手法は連続体モデルに基づく解析であるので，合金に限らず，セラミクスやポリマーアロイの相分離の
図4

A 40 at  B 合金の 1000 K 等温時効における 2次元相分離シミュレーション．
＝0, (b) t′＝10, (c) t′＝20, (d) t′＝50,(a ) t ′(e ) t ′＝100, (f) t′＝200, (g) t′＝300, (h) t′＝500.

ま て り あ 第48巻 第10号(2009)Materia Japan

計算にも適用できる利点がある．特に近年，種々の分野で活発に研究が進められているナノ・メゾスケールにおける材料
設計では，デバイスのスケール自体が微視的な材料組織のス
ケールに重なってくるので，組織形態の解析が，これまでに

も増して重要となりつつある．この時，本稿のようなプログ
ラミングのスキルを身につけておくことは，次世代の新材
料・デバイス設計において，有力な武器になるのではないだ
ろうか．
なお本稿では，より基礎的な事項やさらに進んだ計算理論( 6 ) J. E. Hilliard: Phase Transformation, ed. by H. I. Aaronson,
ASM, Metals Park, Ohio, (1970), 497.
( 7 ) N. Saunders and A. P. Miodownik: CALPHAD, Pergamon,(1998).
( 8 ) J. W. Cahn: The Selected Works of J. W. Cahn, ed. by W. C.
Carter and W. C. Johnaon, TMS, (1998), 29.
( 9 ) 小山敏幸ふぇらむ，11(2006), 647.

の詳細については紙面の関係上ふれなかった．計算理論に関する基礎に関しては文献(1)(3)等をご参照願いたい．また★★★★★★★★★★★★★★★★★★★★★

各種のデモプログラムやより進んだ計算理論の詳細について
は ， 著 者 の ホ ー ム ペ ー ジ ( http: // www.nims.go.jp / mpsg /PhaseField_Modeling.htm)にて公開しているので，興味のある方は参照していただきたい．次回は，マルテンサイト変態の 2 次元シミュレーションプログラムについて解説する
予定である．
小山敏幸

文献( 1 ) 小山敏幸ふぇらむ，9(2004), 240, 301, 376, 497, 905.
( 2 ) 小山敏幸まてりあ，42(2003), 397, 470.
( 3 ) T. Koyama: Chapter 21 in Springer Handbook of Materials
Measurement Methods, H. Czichos, T. Saito and L. Smith(Eds), SpringerVerlag, (2006), 10311054.
( 4 ) T. Koyama: Science and Technology of Advanced Materials,
9(2008), 013006.
( 5 ) T. Miyazaki, A. Takeuchi, T. Koyama and T. Kozakai: Trans.
JIM, 32(1991), 915920.



1988年 3 月 名古屋工業大学大学院 工学研究科 博
士前期課程修了
1990年 4 月 名古屋工業大学大学院 工学研究科 博
士後期課程単位取得後退学
1990年 4 月 名古屋工業大学 工学部 材料工学科
助手
独 物質・材料研究機構 計算材料科学研
2002年 4 月 究センター 主任研究員
独 物質・材料研究機構 計算材料科学研
2005年 4 月 究センター 主幹研究員
独 物質・材料研究機構 計算科学センタ
2006年 4 月 ー 主幹研究員
独 物質・材料研究機構 新構造材料セ2009 年 4 月 
ンター 主幹研究員
現在に至る．
博士(工学)

専門分野材料工学
主に材料組織形成の計算機シミュレーション研究に従事．

★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★

入 門 講 座