はじめに
フェーズフィールド法は材料組織形態の時間変化を予測するシミュレーション手法です。今回はフェーズフィールド法で用いられる最もベーシックな方程式であるアレン-カーン方程式の導出についてまとめます。
アレン-カーン方程式
以下にアレン-カーン方程式を示します。
\frac{\partial \phi}{\partial t}=-M_{\phi}\frac{\delta G}{\delta \phi} \tag{1}\\
\end{align}
ここで\(G\)は時間\(t\)に関する関数で系全体のギブスの自由エネルギーを表します。また、\(\phi\)は秩序変数、\(M_{\phi}\)はフェーズフィールドモビリティと呼ばれる界面の移動しやすさに関するパラメーターです。
アレン-カーン方程式は非保存量の秩序変数\(\phi\)（凝固、相変態、再結晶など）の時間発展を記述する方程式です。（対して、保存量の秩序変数はカーン-ヒリアード方程式で記述されます。）これを解くには材料組織中のギブスの自由エネルギー\(G\)をモデル化し、汎関数微分項\(\frac{\delta G}{\delta \phi}\)を計算する必要があります。
物理的背景と導出
アレン-カーン方程式\((1)\)はギブスの自由エネルギー最小化の概念から導くことが出来ます。ギブスの自由エネルギーの最小化とは「エネルギーの出入りがない閉じた系において、等温・等圧過程の自発変化はギブスの自由エネルギーGが減少する方向に生じる。」という熱力学上の性質です。このことは数式を用いると次のように表現できます。
\frac{dG}{dt}\leq0 \tag{2}\\
\end{align}
ここで式\((2)\)左辺は汎関数微分\(\frac{\delta G}{\delta \phi}\)を用いて、次のように書けます。
\frac{dG}{dt}=\int_{V} \frac{\delta G}{\delta \phi}\frac{\partial \phi}{\partial t}dV\tag{3}\\
\end{align}
唐突に汎関数微分を用いた変形を行ってしまいしたが、これについて詳細を説明しだすと本筋からずれてしまうので、この部分は後から考えることとします。上式を式\((2)\)に代入すると、
\int_{V} \frac{\delta G}{\delta \phi}\frac{\partial \phi}{\partial t}dV\leq0 \tag{4}\\
\end{align}
さて、上式が常に成り立つ場合、秩序変数の時間発展\(\frac{d\phi}{dt}\)はどのような形を取るのが適切でしょうか？もっとも簡単には、以下の関係式がよさそうです。
\frac{\partial \phi}{\partial t}= -M_{\phi}\frac{\delta G}{\delta \phi}\tag{5}\\
\end{align}
ここで\(M_{\phi}\)は正の定数とします。確認のため、式\((5)\)を式\((4)\)に代入してみます。
-M_{\phi}\int_{V} \left(\frac{\delta G}{\delta \phi}\right)^2dV\leq0 \tag{6}\\
\end{align}
左辺の積分は中身が二乗なので符号は正になります。また\(M_{\phi}\)も正の定数なので、これにマイナスがかけられて、左辺は常に０以下になることが分かります。
式\((5)\)はアレン-カーン方程式\((1)\)そのものです。以上より、アレン-カーン方程式はギブスの自由エネルギーを減らしていく方向に、秩序変数\(\phi\)を発展させていく方程式であることが示せました。
汎関数と汎関数微分
さて、前節では汎関数微分\(\frac{\delta G}{\delta \phi}\)を用いた変形を何の説明もなく使ってしまいましたが、そもそも汎関数や汎関数微分ってなんだっけ？という感じです。ここからはこの辺りについて、もろもろおさらいしておきます。
例えば、\(I=F(x)\)という通常の関数があるとしましょう。当然ですが、この時、\(x\)の値が決まると\(I\)の値が決まります。では次のような式はどうでしょうか？
I =\int _{a}^{b}F(x)dx\tag{7}\\
\end{align}
この場合は、右辺は\(x\)に関する定積分なので、\(I\)は\(x\)の関数ではありません。一方で、関数\(F(x)\)の形を変えると、\(I\)の値が変わるので、見方によっては\(I\)は関数\(F(x)\)の関数といってもよいでしょう。このように、ある関数\(F(x)\)の形が決まることで、その値\(I\)が定まるとき、\(I\)は\(F(x)\)の汎関数と呼ばれます。
ここで、\(F(x)\)がさらに別の関数\(f(x)\)の関数である、すなわち\(F(f(x))\)のような形で表される場合、式\((1)\)は次のように書けます。
I =\int _{a}^{b}F(f(x))dx\tag{8}\\
\end{align}
ややこしいですが、状況は特に変わりません。やはり\(I\)は\(F\)の汎関数です。さらにさらに、\(F\)が\(f(x)\)の一階微分\(\frac{df}{dx}\)の関数\(F\left(f,\frac{df}{dx}\right)\)とすれば、
I =\int _{a}^{b}F\left(f,\frac{df}{dx}\right)dx\tag{9}\\
\end{align}
見やすくするため\(x\)の表記は省略しました。だからなんなの？という気がしてきますが、もう少し続きます。今、関数\(F\)の形を少しだけ変えて、汎関数\(I\)の変分\(\delta I\)を求めることを考えます。通常の微分において微小変動を表す際には、\(d\)を用いて\(dI\)などと表しますが、汎関数の場合は「変分」として\(\delta\)を用います。
ここでいう変分\(\delta \)とは\(x\)を\(x\rightarrow x+dx\)と動かすのではなく、あくまでも関数\(F(x)\)の形を少しだけ変えるをことを意味します。例えば、\(F(x)=x^2\)とすると\(F(x)=\lim_{\epsilon \to 0}x^2+\epsilon x\)のようなイメージです。（というかそもそも\(x\)は積分によって値が決まるので、動かすこと自体出来ない。）
逆に言うと式\((7)\sim(9)\)のような関数において、極値探索のために\(I\)の微小変化量を求めようとした場合、\(x\)の微分からそれを見つけることは出来ないため、上記のような「関数形自体を少しだけ変化させる」という考え方が出てくるわけです。
さて、天下りではありますが、上記、汎関数\(I\)の変分\(\delta I\)は最終的に次式にて与えられます。導出は長いのでトグルに記載しておきます。
\delta I &=\int _{a}^{b}\left\{ \frac{\partial F}{\partial f}- \frac{d}{dx}\frac{\partial F}{\partial\left(\frac{df}{dx}\right)}\right\}\delta f dx\\
\tag{10}
\end{align}
式\((9)\)の両辺で変分をとると次式が得られます。
\delta I &=\int _{a}^{b}\delta F\left(f,\frac{df}{dx}\right)dx\\
\tag{A1}
\end{align}
ここで、
F(f(x))の場合、\delta F &=\frac{\partial F}{\partial f}\delta f\tag{A2}
\end{align}
の性質を利用すると、式\((A1)\)は次式のように書けます。
=\int _{a}^{b}\left\{\frac{\partial F}{\partial f}\delta f+ \frac{\partial F}{\partial \left(\frac{df}{dx}\right)}\delta \left(\frac{df}{dx}\right)\right\}dx\tag{A3}
\end{align}
式\((A2)\)では急に偏微分記号が出てきて一瞬混乱しますが、具体例で考えるとわかりやすいです。例えば、\(F=f^2\)という関数を考えると、
\delta F &=(f+\delta f)^2-f^2\\
&=2f\delta f-\delta f^2
\end{align}
変分の二乗項は微小として無視すると、
&=2f\delta f=\frac{\partial F}{\partial f}\delta f
\end{align}
以上のように\(\delta F\)は\(f\)による偏微分を用いた連鎖律で表現可能であることが分かります。
さらに、変分\(\delta\)と微分\(\frac{d}{dx}\)を入れ替えられる性質
\delta\left(\frac{df}{dx}\right) &=\frac{d}{dx}\left(\delta f\right)\tag{A4}
\end{align}
を利用すると、式\((A3)\)は次式のように書けます。
=\int _{a}^{b}\left\{\frac{\partial F}{\partial f}\delta f+ \frac{\partial F}{\partial \left(\frac{df}{dx}\right)}\frac{d \left(\delta f\right)}{dx}\right\}dx\tag{A5}
\end{align}
式\((A4)\)に関しても具体例で考えるとわかりやすいです。例えば、\(f(x)=x^2\)という関数を考えると、
式(A4)左辺&=\delta\left(\frac{df}{dx}\right) =\delta(2x)=2\delta x\\
式(A4)右辺&=\frac{d}{dx}\left(\delta f\right) \\
&=\frac{d}{dx}\left(f(x+\delta x)-f(x )\right)\\
&=\frac{d}{dx}\left\{(x+\delta x)^2-x^2\right\}\\
&=\frac{d}{dx}\left\{2x\delta x+\left(\delta x\right)^2\right\}\\
&\sim \frac{d}{dx}\left(2x\delta x\right)\\
&= 2\delta x+ 2x\frac{d(\delta x)}{dx}\\
&\sim 2\delta x
\end{align}
左辺＝右辺より、式\((A4)\)が成り立つことが示せました。
さらに式\((A5)\)右辺第二項は部分積分を用いて以下のように変形できます。
式(A5)右辺第二項&=\int _{a}^{b}\frac{\partial F}{\partial \left(\frac{df}{dx}\right)}\frac{d \left(\delta f\right)}{dx}dx\\
&=\left[\frac{\partial F}{\partial \left(\frac{df}{dx}\right)} \delta f\right]_{a}^{b}-\int _{a}^{b}\frac{d}{dx}\left\{\frac{\partial F}{\partial \left(\frac{df}{dx}\right)} \right\}\delta fdx\\
&=-\int _{a}^{b}\frac{d}{dx}\left\{\frac{\partial F}{\partial \left(\frac{df}{dx}\right)} \right\}\delta fdx
\tag{A6}
\end{align}
ここで二行目から三行目への変形では、領域境界で変分が０になる仮定、\(\delta f(a)=\delta f(b)=0\)を用いて右辺第一項を消去ました。この仮定は境界条件が定められている状況を意味していますが、これが物理的に妥当かはその問題ごとに考える必要があります。
以上、式\((A6)\)を式\((A5)\)に代入すると
\delta I&=\int _{a}^{b}\left\{\frac{\partial F}{\partial f}\delta f+ \frac{\partial F}{\partial \left(\frac{df}{dx}\right)}\frac{d \left(\delta f\right)}{dx}\right\}dx\\
&=\int _{a}^{b}\left[\frac{\partial F}{\partial f}\delta f-\frac{d}{dx}\left\{\frac{\partial F}{\partial \left(\frac{df}{dx}\right)} \right\}\delta f\right]dx\\
&=\int _{a}^{b}\left\{\frac{\partial F}{\partial f}-\frac{d}{dx}\frac{\partial F}{\partial \left(\frac{df}{dx}\right)} \right\}\delta fdx
\tag{A7}
\end{align}
上式は式\((10)\)と一致しています。
ここで右辺括弧内を汎関数微分\(\frac{\delta I}{\delta f}\)として定義します。
\frac{\delta I}{\delta f} &= \frac{\partial F}{\partial f}- \frac{d}{dx}\frac{\partial F}{\partial\left(\frac{df}{dx}\right)}\tag{11}
\end{align}
この表記を用いると式\((10)\)は次のように書けます。
\delta I &=\int _{a}^{b}\frac{\delta I}{\delta f}\delta f dx
\tag{12}
\end{align}
ここで、\(\delta f\)は\(x\)の関数であり、積分の外には出せないことに注意しましょう。上式は式\((9)\)の全微分形式になっています。すなわち、汎関数微分\(\frac{\delta I}{\delta f}\)は\(f\)の関数形を少し変化させた場合に\(I\)の関数形がどの程度変わるかの変化率を表します。
また、本来、\(F\)は\(f(x)\)以外にも\(f^{\prime}(x)\)を変数を持つにもかかわらず、式\((12)\)に含まれる微小変動量は\(\delta f\)のみであることに注意します。導出過程を確認すると分かりますが、\(\delta f\)を微分したものが\(\delta f^{\prime}\)なので、これらは結果的に\(\delta f\)を用いて表現することができるためです。
アレン-カーン方程式における汎関数微分
実際にアレン-カーン方程式\((1)\)を解くには、ギブスの自由エネルギーの汎関数微分\(\frac{\delta G}{\delta \phi} \)を計算する必要があります。ここまでは汎関数の一般論を述べてきましたが、ここからはこれらの議論がアレンカーン方程式の場合どのように適用されるのか確認します。系全体のギブスの自由エネルギー\(G\)は空間密度\(g\)を用いると次のように書けます。
G(t)=\int_V g(\phi(\boldsymbol{x},t), \nabla\phi(\boldsymbol{x},t))dV\tag{13}
\end{align}
ここで：
\nabla \phi&=\left(\frac{\partial \phi}{\partial x},\frac{\partial\phi}{\partial y},\frac{\partial \phi}{\partial z}\right)\tag{14}\\
\end{align}
上記の通り、\(g\)は秩序変数\(\phi, \nabla\phi\)の関数\(g(\phi, \nabla\phi)\)として書けます。このことはギブスの自由エネルギー密度が相界面の勾配、相界面の曲率に依存することを示しています。（別記事参照）また、\(\phi\)は三次元空間\(\boldsymbol{x}\)および時間\(t\)の関数\(\phi(\boldsymbol{x},t)\)です。右辺は空間積分を含むので、\(G\)は最終的に時間\(t\)のみの関数となります。さて、ここで式\((9)\)を再掲します。
I =\int _{a}^{b}F\left(f,\frac{df}{dx}\right)dx\tag{9：再掲}\\
\end{align}
式\((9)\)は関数\(F\)に関する汎関数\(I\)を表していました。これと式\((13)\)を比較すると、\(G\)は\(g\)に関する汎関数になっていることが分かります。したがって、汎関数\(G\)の変分\(\delta G\)は式\((11)\)とのアナロジーから次のように書けます。
\frac{\delta G}{\delta \phi}&=\frac{\partial g}{\partial \phi}
-\nabla \cdot
\left(\frac{\partial g}{\partial \left(\nabla \phi\right)}\right)
\tag{15}
\end{align}
式\((13)\)の両辺で変分をとると次式が得られます。
\delta G(\phi, \nabla\phi)&=\int_V \delta g(\phi, \nabla\phi)dV\\
&=\int_V \frac{\partial g}{\partial \phi}\delta \phi
+\frac{\partial g}{\partial (\nabla \phi)}\cdot \delta(\nabla\phi)dV
\tag{B1}
\end{align}
ここで右辺第二項はベクトル関数同士の内積になっていることに注意しましょう。この項はベクトル関数による微分を含んでおり少々分かりづらいですが、成分で書き下すと次のように書けます。
\int_V\frac{\partial g}{\partial (\nabla \phi)}\cdot \delta(\nabla\phi)dV
&=
\int_V
\begin{pmatrix}
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial x}\right)}\\
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial y}\right)}\\
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial z}\right)}
\end{pmatrix}
\cdot
\begin{pmatrix}
\delta \left(\frac{\partial \phi}{\partial x}\right)\\
\delta \left(\frac{\partial \phi}{\partial y}\right)\\
\delta \left(\frac{\partial \phi}{\partial z}\right)\\
\end{pmatrix}
dV\\
&= \int_V
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial x}\right)}
\delta \left(\frac{\partial \phi}{\partial x}\right)
+
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial y}\right)}
\delta \left(\frac{\partial \phi}{\partial y}\right)
+
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial z}\right)}
\delta \left(\frac{\partial \phi}{\partial z}\right)dV
\tag{B2}
\end{align}
上式に部分積分を適用すると次のよう変形できます。
&=\int_V
\frac{\partial }{\partial x}
\left(
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial x}\right)}
\delta \phi
\right)
–
\frac{\partial }{\partial x}
\left(
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial x}\right)}
\right)
\delta \phi \\
&\qquad \qquad
+
\frac{\partial }{\partial y}
\left(
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial y}\right)}
\delta \phi
\right)
–
\frac{\partial }{\partial y}
\left(
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial y}\right)}
\right)
\delta \phi\\
&\qquad \qquad \qquad+
\frac{\partial }{\partial y}
\left(
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial z}\right)}
\delta \phi
\right)
–
\frac{\partial }{\partial z}
\left(
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial z}\right)}
\right)
\delta \phi dV\\
&= \int_V
\nabla \cdot
\begin{pmatrix}
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial x}\right)}\delta \phi\\
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial y}\right)}\delta \phi\\
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial z}\right)}\delta \phi
\end{pmatrix}
-\nabla \cdot
\begin{pmatrix}
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial x}\right)}\\
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial y}\right)}\\
\frac{\partial g}{\partial \left(\frac{\partial \phi}{\partial z}\right)}
\end{pmatrix}
\delta \phi
dV\\
&= \int_V
\nabla \cdot
\left(\frac{\partial g}{\partial \left(\nabla \phi\right)}\delta \phi\right)
-\nabla \cdot
\left(\frac{\partial g}{\partial \left(\nabla \phi\right)}\right)
\delta \phi
dV
\tag{B3}
\end{align}
第一項についてガウスの発散定理を用いると次のように書けます。
= \int_S
\left(\frac{\partial g}{\partial \left(\nabla \phi\right)}\delta \phi\right)\cdot \boldsymbol{n} dS
-\int_V \nabla \cdot
\left(\frac{\partial g}{\partial \left(\nabla \phi\right)}\right)
\delta \phi
dV\tag{B4}\\
\end{align}
上記の\(\boldsymbol{n}\)は境界に垂直な方向（法線方向）の単位ベクトルです。今、領域境界では秩序変数\(\phi\)の値が固定されている（\(\phi = const\)、ディリクレ条件に相当。）と仮定すると、\(\delta \phi = 0\)と式\((B4)\)の第一項は0となり、最終的に式\((B1)\)第二項は次のように書けます。
=-\int_V \nabla \cdot
\left(\frac{\partial g}{\partial \left(\nabla \phi\right)}\right)
\delta \phi dV\tag{B5}\\
\end{align}
今回、式\((B4)\)の第一項を消去するために、\(\delta \phi = 0\)の仮定を用いましたが、実際に問題を解く際に、このような仮定が成立する状況は多くはありません。どちらかというと\(\nabla \phi \cdot \boldsymbol{n}=0\)(秩序変数の出入りがない閉じた系、ノイマン条件に相当。)を仮定し、式\((B4)\)第一項を０とする場合が多いようです。なぜなら\(g\)の具体的な形として、勾配エネルギー密度\(\frac{a^2}{2}(\nabla \phi)^2 \)を適用すると、
\int_S
\left(\frac{\partial g}{\partial \left(\nabla \phi\right)}\delta \phi\right)\cdot \boldsymbol{n} dS=\int_S a^2\delta \phi \left(\nabla\phi \cdot \boldsymbol{n}\right) dS
\end{align}
が成り立ち、\(\nabla \phi \cdot \boldsymbol{n}=0\)の条件でも式\((B4)\)第一項を消去できるためです。（勾配エネルギー密度は次回記事で導入します。）
以上、式\((B5)\)を式\((B1)\)に代入すると
\delta G(\phi, \nabla\phi)&=\int_V \delta g(\phi, \nabla\phi)dV\\
&=\int_V \frac{\partial g}{\partial \phi}\delta \phi
+\frac{\partial g}{\partial (\nabla \phi)}\cdot \delta(\nabla\phi)dV\\
&=\int_V \frac{\partial g}{\partial \phi}\delta \phi
-\nabla \cdot
\left(\frac{\partial g}{\partial \left(\nabla \phi\right)}\right)
\delta \phi dV\\
&=\int_V \left\{\frac{\partial g}{\partial \phi}
-\nabla \cdot
\left(\frac{\partial g}{\partial \left(\nabla \phi\right)}\right)
\right\}
\delta \phi dV
\tag{B6}
\end{align}
上式について、汎関数微分の定義、
\delta G(\phi, \nabla\phi)=\int_V \frac{\delta G}{\delta \phi}
\delta \phi dV
\tag{B7}
\end{align}
と比較すると、
\frac{\delta G}{\delta \phi}=\frac{\partial g}{\partial \phi}
-\nabla \cdot
\left(\frac{\partial g}{\partial \left(\nabla \phi\right)}\right)
\tag{B8}
\end{align}
となり、式\((15)\)を得ることが出来ました。
以上でギブスの自由エネルギーの汎関数微分\(\frac{\delta G}{\delta \phi} \)が導出できました。ただし実際にこれを計算する場合には、空間エネルギー密度\(g\)の形を定義する必要があります。これについては、次回の記事でまとめたいと思います。
おわりに
今回はアレン-カーン方程式の導出についてまとめました。フェーズフィールド法において汎関数の概念は非常に重要ですが、直感的に理解しづらい部分も多いように思います。特に式変形時に適用される仮定は、教科書だとさらっと導入されますが、一つ一つかみ砕いていかないと理解できないと感じました。また、アレン-カーン方程式自体は導出もシンプルで非常に汎用的な式ですが、それゆえに自由エネルギー密度\(g\)を如何にモデル化するかが予測精度を左右すると考えられます。次回はこの\(g\)の計算方法についてまとめたいと思います。
参考文献
[1] 高木 知弘ら, “フェーズフィールド法”, 養賢堂, 2012/3/2
[2] 山中 晃徳ら, “Pythonによるフェーズフィールド法入門 基礎理論からデータ同化の実装まで”, 丸善出版, 2023/12/15
コメント