熱物性論 2015 （松本充弘）
： p. 33
非平衡系を扱う：TDGL モデル
6
前章では，相転移がエネルギーとエントロピーの競合によって起きることを述べた．そ
の概略を，グラフで復習すると下図のようになる．我々はこの 平均場近似モデル に基づ
いて，平衡状態を議論した．すなわち，完全な平衡状態では，ミクロ自由エネルギー f が
極小となる状態（上図の ×）が実現すると考えた．この章では × 以外の状態にある系，つ
まり非平衡状態の系，をどう扱えばいいかについて考えよう．
Entropy S
Energy E
Free Energy
F = E - TS
T < Tc
X
X
T > Tc
Compositon
Compositon
X
Compositon
図 6–10: エネルギーとエントロピーの競合による相転移のメカニズムの説明．
例えば，高温で均一に混ざっている A-B ２元合金系を，相転移温度以下に冷やすと，A-rich
な領域と B-rich な領域に次第に分離していくことを，モンテカルロ計算で眺めることがで
きる．このような 非平衡状態からの相変化の進行 に伴って，特徴的な空間構造が生まれ
ることが多い．本章の目的は，相変化に伴う構造形成の過程 (構造の時間変化) を，平均場
描像に基づく簡単な自由エネルギーモデルによって調べることである．
6.1
自由エネルギーのモデル
相変化を議論するのに必要最小限の性質をもつ簡単なモデルを，前章の平均場近似の結
果を参考にして構築しよう．
• Step 1: 状態を表すパラメタ（秩序変数，order parameter）s を定義する．一般に， 秩序変数としては，例えば，合金系な
s は場所ごとに決まるので，空間座標 ~r の関数 s(~r) である．習慣的に，s = 0 が完全
ら組成，気液相転移なら密度，磁性転
移なら磁化の強さなどが選ばれる．
無秩序状態を表すように定義することが多い．例えば，合金系なら，組成 c に対し
て，s = 2c − 1 とすればよい．
• Step 2: 自由エネルギー密度 f (s) を，s のべき関数で近似する．最も簡単なモデル
では，左右対称（偶関数）でよいので，奇数べきの項は必要ない．また，エネルギー
の原点は自由に選べるので，定数項はゼロとする．したがって
f (s) = Bs2 + Ds4 + F s6 + · · ·
ここで，係数 B, D, F, . . . は温度に依存してもよい．
• Step 3: 低温において，平均場近似モデルのように３つの極値 (２つの極小とその
間の極大) をもつためには，少なくとも４次関数である必要がある．
• Step 4: 相転移温度 T = Tc で３つの極値が１つに融合するためには，２次の係数
に温度依存性をもたせて，B(T = Tc ) = 0 とすればよい．
平均場近似では，エントロピー項が
s log s の形で表されていた．このまま
でも悪くはないが，モデルとしてはべ
き関数だけで表されている方が計算が
簡単である．また，対数関数は s → 0
などで発散するから，取り扱いがやや
面倒でもあり，単純なべき関数のほう
が好まれる．
熱物性論 2015 （松本充弘）
： p. 34
これらの条件を満たす最も簡単な自由エネルギーモデルが，次式である：
f (s) = b · (T − Tc ) s2 + D s4
(57)
ここで，b と D は正の定数と考える．
系全体の自由エネルギーは，この局所的な自由エネルギー密度 f (s) を空間積分すること
によって得られると考える：
∫
F [s] =
[
]
b · (T − Tc ) s(~r)2 + D s(~r)4 dV
(58)
すなわち，F は s(~r) の汎関数 functional と考える．
6.2
界面自由エネルギー
相転移温度より低温で，２つの相が分離・共存しているとき，その界面は余分な自由エ このため，界面は可能な限り，常に小
さくなろうとする．これが表面張力
ネルギーを持っていると考えられる．「界面」とは，s が大きく変化しているところであ surface tension の由来である．後の
~ 2 に比例する 章でもう少し詳しく考察する予定で
るから，界面が持つ余分な自由エネルギーを，s の勾配の２乗，つまり |∇s|
ある．
項としてモデルに取り入れる：
∫ [
F [s(~r)] =
]
~ 2 dV
b · (T − Tc ) s(~r)2 + D s(~r)4 + σ|∇s|
(59)
ここで，σ は表面張力の大きさに対応する正定数である．
勾配を考慮しなければ，界面は数学
的にいくらでも薄くなる（このとき，
∇s → ∞）ことが可能になる．これ
では，現実に観察される「有限の厚さ
を持った界面」とは全く異なった非現
実的なモデルになってしまう．
この形の自由エネルギーモデルを，Ginzburg-Landau 型自由エネルギー といい，もとも
とは超伝導機構の現象論的説明のために考えられた有名なモデルであり，相転移ダイナミ Vitaly L. Ginzburg (1916–2009) ロ
クスを議論する出発点として今でもよく用いられている．また，この自由エネルギーをさ シアの理論物理学者．超伝導・超流動
の研究により，2003 年ノーベル物理
まざまな形で拡張したものは，混相流や材料学の数値計算においてフェーズフィールドモ 学賞．
デル phase field model として広く用いられている．
6.2.1
Lev Davidovich Landau (1908–1968)
ロシアの理論物理学者．ヘリウム超流
動の研究により，1962 年ノーベル物理
学賞．その業績は多岐に渡る．Evgeny
Mikhailovich Lifshitz らと共同で執
筆した教科書「理論物理学教程」は有
名．
(補足) Ginzburg-Landau モデル と Cahn-Hilliard モデル
Ginzburg と Landau が，超伝導相転移を説明するマクロスコピック理論をロシア語の
雑誌（Zh. Eksp. Teor. Fiz.）に発表したのは，1950 年である．おそらくそれとは独立
に (当時の冷戦下では，科学の分野においてもソ連・共産圏と西欧諸国の交流は難しかっ
た)，非一様な溶液系の自由エネルギーモデルを，Cahn と Hilliard が 1958 年に提案して
いる [J.W. Cahn, J. Hilliard, “Free Energy of a Nonuniform System. I. Interfacial Free
Energy” J. Chem. Phys., 58 (1958) 258–267]．
∫ [
]
2
F =
f0 + κ (∇c) + · · · dV
(Wikipedia)
(60)
c は溶液の組成である．界面の効果を組成勾配の２乗で表現するという意味において，これ
は明らかに，式 (59) と本質的に同じである．このような歴史的背景から，Ginzburg-Landau
モデル（主に物性物理学や素粒子物理学の分野が多い気がする），Cahn-Hilliard モデル
（化学の分野でよく見かけるほか，流体物理学などでも使われる） の両方の呼称が使われ
ている．
熱物性論 2015 （松本充弘）
： p. 35
6.2.2
(参考) Phase field 法
Ginzburg-Landau 型の自由エネルギーモデルから出発して，相界面の時間発展を記述するモデル
が phase field model の名前で提案され，材料工学や混相流計算の分野で広く使われるようになっ
た [A.A. Wheeler, W.J. Boettinger, G.B. McFadden, “Phase-field model for isothermal phase
transitions in binary alloys,” Phys. Rev. A 45 (1992) 7424–7440 および R. Kobayashi, “Modeling
and numerical simulations of dendritic crystal growth,” Physica D 63 (1993) 410–423]. 以下， http://tkoyama.web.nitech.ac.jp/docs
/Lecture H20/H20 Chapter 11.pdf
小山敏幸氏の解説より，かいつまんで紹介する：
まず，phase-field 変数とよばれるスカラー連続変数 φ (~
r) (0 ≤ φ ≤ 1) を局所的に定義する．こ
れは，対象とする２つの相を表現するもので，例えば「合金融液からの凝固過程」を扱う場合には
{
φ=
∼0
∼1
液相
固相
(61)
のように定義すればよい．系全体の自由エネルギー F を次の形で定義する：
∫ [
F [φ] =
]
f (φ, c, T ) +
1 2
|∇φ|2 dV
2
(62)
ここで，組成 c に依存するバルク自由エネルギー密度は，例えば
f (φ, c, T ) = h(φ)f S (c, t) + [1 − h(φ)]f L (c, T ) + W g(φ)
(63)
のようにモデル化する．これは，固相と液相の自由エネルギー密度 f S と f L の線形和に，中間濃度
の「ハンディキャップ」W を加えたモデルであり，その重み関数としては，例えば
h(φ)
g(φ)
(
=
φ3 10 − 15φ + 6φ2
=
φ − 2 (1 − φ)
)
2
(64)
(65)
といった形がよく用いられるそうである．明らかに，式 (62) は Ginzburg-Landau 型の自由エネル
ギーである．
凝固問題（結晶化過程）を扱うときの特徴は，界面エネルギーの係数 に 結晶方位依存性 をもた
せることである．たとえば
(θ) = 0 [1 + ν cos(kθ)]
(66)
によって，方位 θ への依存性がモデル化されている．この解説 [右欄注釈] には，Ode らによる，
Fe-0.5 mol%C 合金のデンドライト成長の計算例 [M. Ode and T. Suzuki, ISIJ International, 42
(2002) 368] が掲載されている．
図 6–11: Phase field 法によって計算した Fe-C 合金の初期凝固デンドライト成長 [Ode and
Suzuki, 2002]. (a) 1.0 ms, (b) 3.0ms, (c) 5.0ms, (d) 10ms, (e) 17ms, (f) 24ms.
熱物性論 2015 （松本充弘）
： p. 36
(参考) 混相流への phase-field kodel の応用例
高田尚樹氏@産総研のホームページから
Phase field model による混相流の
計算例：液滴の自由落下と液膜への
衝突．同じく高田尚樹氏@産総研の
ホームページから
↓
↓
↓
↓
熱物性論 2015 （松本充弘）
： p. 37
最近見つけたネットワーク理論に関する面白い本から．2005 年出版．ランダウはいかに天才的な仕事をしたことか…
熱物性論 2015 （松本充弘）
： p. 38
途中省略
熱物性論 2015 （松本充弘）
： p. 39
熱物性論 2015 （松本充弘）
： p. 40
6.3
非平衡状態からの緩和 – TDGL 方程式
系が，f (s) の極小点以外の状態（＝非平衡状態）にある時は，時間がたつにつれて極小
点に向かって緩和していくと考えられる．直感的（経験的）に，その緩和の速さは 一般化
df
力 −
に比例すると考えられるので，各場所において，次のような運動方程式（時間発
ds
展方程式 time-development eq.）モデルを仮定する：
ds
δF
δ
= −k
（ は汎関数微分）
dt
δs
δs
(67)
ここで k は s の変わりにくさ（「慣性」inertia ）をあらわす適当な正の定数である．F の
モデル，(59) 式，について汎関数微分を実行すると
ds
= −2kb(T − Tc )s − 4kDs3 + 2kσ∆s
dt
勾配の２乗の項の汎関数微分は，
∫
(68)
∫
~ + δs)|2 dV −
|∇(s
∫
=2
という s(~r) の時間発展方程式が得られる．さらに，時間，空間，温度の単位を適当に無次
元化 すると
~ 2 dV
|∇s|
~ ∇δsdV
~
∇s
+ o(δs)
を部分積分することによって Lapla-
ds
= −T ∗ s − s3 + ξ∆s
dt
cian ∆ ≡
(69)
∂2
∂~
r2
で表現できる．
式 (68) には，k, b, D, σ の４つのパラ
が得られる．この式は，TDGL (Time-dependent Ginzburg-Landau) 方程式とよばれ，非 メタが含まれている．一方，時間・空
間・温度の３つの単位があるから，適
平衡系の構造形成を議論する出発点としてよく用いられるモデルである．ここで，
T ∗ ∝ T − Tc
(70)
当な単位系を使って無次元化すること
で，パラメタを１つに減らすことがで
きる．式 (69) のように，表面張力に相
当するパラメタ ξ を残すことが多い．
は，相転移温度を基準として測った無次元化温度であり，相転移温度からのずれを表すの
で，負の値もとり得ることに注意してほしい．また，ξ は無次元化された表面張力パラメ
タである．
6.4
TDGL 方程式の数値計算
式 (69) は，ポテンシャル
f (s) =
T∗ 2 1 4
s + s
2
4
(71)
内の「s の運動方程式」に，界面項が加わったものと見ることができる．
この見方に立って，解の性質をまずは概観しよう：
ポテンシャルの部分は，局所的な s(~
r)
のみで与えられているのに対して，界
面項 ξ∆s は，微分という非局所的な
相互作用である．この非局所性が空間
的な構造形成のカギとなる．
• T ∗ > 0，すなわち相転移温度より高温では，s = 0 が式 (69) の唯一の安定解である．
すなわち，平衡状態で系全体が均一に混ざっている．
• T ∗ < 0 では，s = 0 の解は f の極大に対応する不安定解であり，安定解は
√
± −T ∗
空間的に均一
√
s=
(72)
r − r0
∗
± −T tanh
r0 近傍に界面を持つ
D
√
2ξ
という２つのタイプが存在する．ここで，D ≡
は 界面の厚さ であり，界面
−T ∗
の場所 r0 は任意である．
d2
dx2
tanh x = −2 tanh x + 2 tanh3 x
√
だから，y = tanh 2x は
y − y3 +
d2 y
=0
dx2
の解となり，適当に数係数を合せるこ
とで，式 (69) の時間変化しない一次
元解が得られる．T ∗ < 0 であること
に注意．
熱物性論 2015 （松本充弘）
： p. 41
#define NUM_CELL 100
#define NUM_STEP 5000
#define NUM_SAVE 10
double state[NUM_CELL];
void initial( )
{
int icell;
for (icell=0;icell<NUM_CELL;icell++) {
state[icell]=0.1*(rand( )/(double)RAND_MAX-0.5);
}
}
T= +0.1
0.1
Order Parameter s
#include <stdio.h>
#include <stdlib.h>
0
10
30
100
200
300
500
0.05
0
-0.05
-0.1
0
void data_out(int iout)
{
int icell;
char filen[100];
FILE *filep;
}
int main()
{
int istep,icell,icell0,icell1,iout;
double temp,st,st0,st1;
double dstate[NUM_CELL];
80
100
T= -0.1
0
10
100
500
1000
1500
2000
3000
sqrt(-T)
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
0
20
double xi=1.0;
double dtime=0.05;
scanf("%lf",&temp);
initial();
istep=0; data_out(istep);
40
60
80
Spatial Coordinate
100
T= -1.0
1
Order Parameter s
for (istep=1;istep<=NUM_STEP;istep++) {
for (icell=0;icell<NUM_CELL;icell++) {
icell0=icell-1; if (icell0==-1) icell0=NUM_CELL-1;
icell1=icell+1; if (icell1==NUM_CELL) icell1=0;
st=state[icell];
dstate[icell] = -temp*st
-st*st*st
+xi*(state[icell0]+state[icell1]-2*st);
}
for (icell=0;icell<NUM_CELL;icell++) {
state[icell] += dtime*dstate[icell];
}
if (istep%NUM_SAVE==0) data_out(istep);
}
return 0;
40
60
Spatial Coordinate
0.4
Order Parameter s
sprintf(filen,"%5.5d.dat",iout); printf("%s\n",filen);
filep=fopen(filen,"w");
for (icell=0;icell<NUM_CELL;icell++) {
fprintf(filep,"%7.4f\n",state[icell]);
}
fclose(filep);
20
0
10
30
70
100
500
2000
5000
sqrt(-T)
0.5
0
-0.5
-1
0
20
40
60
80
Spatial Coordinate
100
}
List 6-1: １次元 TDGL 数値計算プログラム例．初期条件として，乱数
図 6–12: １次元 TDGL 数値計算結果の例
を使ったノイズを与えている．
数値計算により時間変化の詳細を追跡できる．ここでは，素朴な差分法 に基づく解法を
とることにする．リスト 6.1 に１次元 TDGL 方程式の数値計算プログラムを，図 6–12
にその結果例を載せた．初期条件として与えた空間的なノイズが，T ∗ > 0 では時間の経過
と共に減衰するのに対し，T ∗ < 0 ではノイズがどんどん成長し，はっきりとした界面が形
成される様子が見える．
熱物性論 2015 （松本充弘）
： p. 42
6.5
ちょっと横道：拡散方程式の数値解法
物質拡散や熱伝導などの輸送現象を表す放物型偏微分方程式 parabolic partial differential
equation (いわゆる 拡散方程式)
以下，記述を簡単にするために１次元
で表記するが，多次元でも同様である．
∂φ(x, t)
∂ 2 φ(x, t)
=D
∂t
∂x2
(73)
の数値解法は，よく知られているように２種類に大別される：
(1) 陽的解法 explicit scheme
右辺の空間差分を時刻 t において行う．すなわち
φ(x, t + ∆t) − φ(x, t)
φ(x + ∆x, t) − 2φ(x, t) + φ(x − ∆x, t)
=D
∆t
(∆x)2
(74)
これを整理して，時刻 t + ∆t での物理量を左辺に集めると
[
]
D∆t
D∆t
D∆t
φ(x, t+∆t) =
φ(x+∆x, t)+ 1 − 2
φ(x, t)+
φ(x−∆x, t) (75)
(∆x)2
(∆x)2
(∆x)2
D∆t
とするときと，次の時刻の φ は，隣り合う格子点で，重み
(∆x)2
r : 1 − 2r : r をつけて φ の空間平均を行うことで求められると解釈できる．１ステッ
無次元量を r ≡
プあたりの計算量が少なく手軽にプログラムを作成できるが，計算の安定性に注意
が必要である．安定に計算が進行するためには，0 < 1 − 2r，つまり
r<
1
2
(76)
を満たす必要がある．List 6.1 で示したプログラムは，この方法を用いている．
(2) 陰的解法 implicit scheme
陽的解法とは異なり，右辺の空間差分を時刻 t + ∆t において行う．
すなわち
右辺の差分を，t と t+∆t を適当に混合
して行うことも可能であり，Crank–
Nicolson 法として知られている．
φ(x, t + ∆t) − φ(x, t)
φ(x + ∆x, t + ∆t) − 2φ(x, t + ∆t) + φ(x − ∆x, t + ∆t)
=D
∆t
(∆x)2
(77)
これを，上と同様に整理すると，
[
]
D∆t
D∆t
D∆t
−
φ(x+∆x,
t+∆t)+
1
+
2
φ(x, t+∆t)−
φ(x−∆x, t+∆t) = φ(x, t)
(∆x)2
(∆x)2
(∆x)2
(78)
従って，次の時刻の φ を求めるためにはこの連立１次方程式を解く必要があり，プ
D∆t
が相当に大きくても安定
ログラミングにも計算にも手間がかかる反面，r ≡
(∆x)2
に計算を進められることが知られている．
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define
#define
#define
#define
NUM_CELL
NUM_STEP
NUM_SAVE
dtime
100
8000
20
0.1
double state[NUM_CELL][NUM_CELL];
FILE *fplot;
void initial() //----------------------------------------{
int ic,jc;
for (ic=0;ic<NUM_CELL;ic++)
for (jc=0;jc<NUM_CELL;jc++)
state[ic][jc]=1.0*(rand()/(double)RAND_MAX-0.5);
}
熱物性論 2015 （松本充弘）
： p. 43
void data_out(int istep) //------------------------------{
int ic,jc;
char cfile[100];
static int flag=0;
FILE *filep;
sprintf(cfile,"%6.6d.dat",istep);
filep=fopen(cfile,"w");
for (ic=0;ic<NUM_CELL;ic++) {
for (jc=0;jc<NUM_CELL;jc++) {
fprintf(filep," %7.4f",state[ic][jc]);
}
fprintf(filep,"\n");
}
fclose(filep);
fprintf(fplot,"splot \’%s\’ matrix\n",cfile);
if (flag==0) {
fprintf(fplot,"pause -1\n");
flag=1;
} else {
fprintf(fplot,"pause 0.01\n");
}
// gnuplot コマンドファイル出力
}
int main() //--------------------------------------------{
int istep,ic,jc,ic0,ic1,jc0,jc1;
double temp,st,grad;
double dstate[NUM_CELL][NUM_CELL];
double xi=1.0;
// 表面張力係数 ξ は 1.0 に固定してある
printf("Input T: ");
scanf("%lf",&temp);
fplot=fopen("tdgl.plt","w"); // gnuplot (Ver 4.6) コマンドファイル作成, pm3d カラーマップを利用
fprintf(fplot,"unset surface\n");
fprintf(fplot,"set style data pm3d\n");
fprintf(fplot,"set view map\n");
fprintf(fplot,"set pm3d implicit at b\n"); // このあたり，古い gnuplot では動かないかも
fprintf(fplot,"set palette rgbformulae 31, -11, 32\n");
fprintf(fplot,"set key at %3d,%3d\n",NUM_CELL,NUM_CELL+10);
fprintf(fplot,"set cbrange [%6.1f:%6.1f] noreverse nowriteback\n",-1.2*fabs(temp),1.2*fabs(temp));
initial();
istep=0;
data_out(istep);
for (istep=1;istep<=NUM_STEP;istep++) {
printf("%5d step\n",istep);
for (ic=0;ic<NUM_CELL;ic++) {
for (jc=0;jc<NUM_CELL;jc++) {
st=state[ic][jc];
ic0=ic-1; if (ic0==-1) ic0=NUM_CELL-1;
ic1=ic+1; if (ic1==NUM_CELL) ic1=0;
grad=state[ic0][jc]+state[ic1][jc]-2*st;
jc0=jc-1; if (jc0==-1) jc0=NUM_CELL-1;
jc1=jc+1; if (jc1==NUM_CELL) jc1=0;
grad+=state[ic][jc0]+state[ic][jc1]-2*st;
dstate[ic][jc]= -temp*st - st*st*st + xi*grad;
}}
for (ic=0;ic<NUM_CELL;ic++) {
for (jc=0;jc<NUM_CELL;jc++) {
state[ic][jc] += dtime*dstate[ic][jc];
}}
if (istep%NUM_SAVE==0) data_out(istep);
}
fclose(fplot);
return 0;
}
List 6.2: ２次元 TDGL 数値計算プログラム例： gnuplot で濃淡図を描くコマンドファ
イルを生成する．古い version では動かないかも知れないので注意．
図 6–13: ２次元 TDGL 方程式の数値計算（List 6.2）の可視化例：T ∗ = −1.0 の場合.
© Copyright 2024 ExpyDoc
