はじめに
前回は以下の界面異方性を考慮したアレン-カーン方程式についてシミュレーションしました。
\frac{\partial \phi}{\partial t}&=-M_{\phi}\left[\left(g_{a}-g_{b}\right)\frac{\partial p(\phi)}{\partial \phi}+W\frac{\partial q(\phi)}{\partial \phi}
-2a\left(\frac{\partial a}{\partial x} \frac{\partial \phi }{\partial x}+\frac{\partial a}{\partial y} \frac{\partial \phi }{\partial y} \right)\right.\\
&\left.\quad\quad\quad\quad\quad\quad-a^2\left(\frac{\partial ^2\phi}{\partial x^2}+\frac{\partial ^2\phi}{\partial y^2} \right)
+\frac{\partial }{\partial x}\left( a\frac{\partial a }{\partial \theta }\frac{\partial \phi }{\partial y}
\right)-
\frac{\partial }{\partial y}
\left(
a\frac{\partial a }{\partial \theta } \frac{\partial \phi }{\partial x}
\right)
\right] \tag{1}\\
a(\theta)&=\bar{a}\left[1+\xi \cos \left\{ k\left( \theta -\theta _0\right)\right\} \right]
\tag{2}\\
q(\phi)&=\phi ^2 (1-\phi)^2 \tag{3}\\
p(\phi)&=\phi ^2(3-2\phi) \tag{4}
\end{align}
ここまでの内容はどちらかというと方程式の性質を理解するためのパラメータースタディが中心でしたが、ここからは実現象を想定したモデル化・パラメーターセットで計算していきます。
その中でも今回は「純物質の凝固」のシミュレーションにフォーカスします。「純物質の凝固」とは、例えばある単体の金属が液体から固体に変化するような状況に相当します。凝固はその点が融点を下回った際に生じる現象なので、温度場も解く必要が出てきます。この時、フェーズフィールドと温度場がどのようにモデル化されるのか、順にまとめていきたいと思います。
エネルギー密度分布関数\(p(\phi)\)の変更
前準備として、エネルギー密度分布関数\(p(\phi)\)を変更します。前回記事では\(p(\phi)\)として式\((4)\)を用いましたが、今回は文献[1]に基づき、以下を用いることとします。
p(\phi)&=\phi ^3(10-15\phi+6\phi^2) \tag{5}\\
\frac{\partial p(\phi)}{\partial \phi}&=30\phi^2(1-\phi)^2\tag{6}
\end{align}
これによりアレン-カーン方程式\((1)\)は次のように修正されます。
\frac{\partial \phi}{\partial t}&=-M_{\phi}\left[30\left(g_{a}-g_{b}\right)\phi^2\left(1-\phi\right)^2+2W\phi (1-\phi)(1-2\phi)
\right.\\
&\left. \quad\quad\quad\quad\quad-2a\left(\frac{\partial a}{\partial x} \frac{\partial \phi }{\partial x}+\frac{\partial a}{\partial y} \frac{\partial \phi }{\partial y} \right)-a^2\left(\frac{\partial ^2\phi}{\partial x^2}+\frac{\partial ^2\phi}{\partial y^2} \right)
\right.\\
&\left.\quad\quad\quad\quad\quad\quad
-\frac{\partial }{\partial x}\left( a\bar{a}k\xi \sin \left\{ k\left( \theta -\theta _0\right)\right\} \frac{\partial \phi }{\partial y}
\right)+
\frac{\partial }{\partial y}
\left(
a\bar{a}k\xi \sin \left\{ k\left( \theta -\theta _0\right)\right\} \frac{\partial \phi }{\partial x}
\right)
\right]
\tag{7}
\end{align}
右辺第一項が前回から変化していることに注意しましょう。
文献[1]ではエネルギー密度分布関数\(p(\phi)\)としてなぜ式\((5)\)を用いているのか、明確な記載はありませんでしたが、エネルギー密度関数の違いは、解に対する化学的駆動力\(\Delta g\)の寄与の仕方に影響を与えます。（過去記事参照）
また式\((4)\)で実際計算してみると、フェーズフィールド場中に不自然な分布が現れてしまうようでした。このあたりの理解はまだ不十分かもしれませんが、ひとまず今回は文献[1]に従っておくこととします。
アレン-カーン方程式への温度場の導入
次にアレンカーン方程式に液相と固相の概念を取り入れます。前回までは\(\phi = 1\)をa相、\(\phi = 0\)をb相としましたが、今回はより具体的に\(\phi = 1\)を固相、\(\phi = 0\)を液相として式\((7)\)を以下のように書き換えます。
\frac{\partial \phi}{\partial t}&=-M_{\phi}\left[30\left(g_{S}-g_{L}\right)\phi^2\left(1-\phi\right)^2+2W\phi (1-\phi)(1-2\phi)
\right.\\
&\left. \quad\quad\quad\quad\quad-2a\left(\frac{\partial a}{\partial x} \frac{\partial \phi }{\partial x}+\frac{\partial a}{\partial y} \frac{\partial \phi }{\partial y} \right)-a^2\left(\frac{\partial ^2\phi}{\partial x^2}+\frac{\partial ^2\phi}{\partial y^2} \right)
\right.\\
&\left.\quad\quad\quad\quad\quad\quad
-\frac{\partial }{\partial x}\left( a\bar{a}k\xi \sin \left\{ k\left( \theta -\theta _0\right)\right\} \frac{\partial \phi }{\partial y}
\right)+
\frac{\partial }{\partial y}
\left(
a\bar{a}k\xi \sin \left\{ k\left( \theta -\theta _0\right)\right\} \frac{\partial \phi }{\partial x}
\right)
\right]
\tag{8}
\end{align}
化学的自由エネルギー密度\(g_{a},g_{b}\)が\(g_{S},g_{L}\)と書き換わっていることに注意しましょう。
凝固が進行する際に、固液界面を移動させる駆動力は化学的自由エネルギー密度の差\(\Delta g =g_{S}-g_{L}\)です。いままで化学的自由エネルギー密度が小さいほうから大きいほうへと界面が移動していく様子をシミュレーションしてきましたが、これが凝固現象そのものとなります。
言い換えると、凝固が進行する際（温度が融点以下の際）、化学的自由エネルギー密度は液相より固相のほうが小さい状態と言えます。溶融が進行する際（温度が融点以上の際）はその逆です。したがって、化学的自由エネルギー密度\(\Delta g \)の符号は温度によって変化します。
では\(\Delta g\)は具体的にどのような形で与えられるのでしょうか？今一度ギブスの自由エネルギーの定義に立ち返りましょう。
液相、固相のギブスの自由エネルギー\(G_{L},G_{S}\rm{[J]}\)は定義より次のように書けます。
G_{L}&=H_{L}-TS_{L}
\tag{9}\\
G_{S}&=H_{S}-TS_{S}\tag{10}\\
\end{align}
ここで、\(H_{L},H_{S}\)は液相、固相のエンタルピー、\(S_{L},S_{S}\)は液相、固相のエントロピーです。
エンタルピーとは等温・等圧過程における熱量そのものを示します。したがって、\(H_{L}\)と\(H_{S}\)は原子をバラバラにするだけのエネルギー、すなわち潜熱\(L\rm{[J]}\)の分だけ異なると考えてよいでしょう。
H_{L}=H_{S}+L
\tag{11}\\
\end{align}
同様に融点\(T_{m}\)における潜熱\(L\)の出入りはエントロピーの変化\(\frac{L}{T_{m}}\)を伴うので、液相のエントロピーは固相のそれよりも\(\frac{L}{T_{m}}\)だけ高いはずです。
S_{L}=S_{S}+\frac{L}{T_{m}}
\tag{12}\\
\end{align}
以上、式\((11),(12)\)を式\((9)\)に代入すると\(G_{L},G_{S}\)の間には次の関係が成り立ちます。
G_{L}&=H_{S}+L-T\left(S_{S}+\frac{L}{T_{m}}\right)\\
&=H_{S}-TS_{S}-L\left(\frac{T}{T_{m}}-1\right)\\
&=G_{S}-L\left(\frac{T-T_{m}}{T_{m}}\right)
\tag{13}\\
\end{align}
この両辺を体積\(V\)で割り、ギブスの自由エネルギー密度\(g_{L},g_{S}\rm{[J/m^3]}\)として考えると
&\frac{G_{L}}{V}=\frac{G_{S}}{V}-\frac{L}{V}\left(\frac{T-T_{m}}{T_{m}}\right)\\
&\rightarrow g_{L}=g_{S}-l\left(\frac{T-T_{m}}{T_{m}}\right)\\
&\rightarrow \Delta g=g_{S}-g_{L}=l\left(\frac{T-T_{m}}{T_{m}}\right)\tag{14}\\
\end{align}
ここで\(l\)は潜熱密度\(\frac{L}{V}\rm{[J/m^3]}\)を表します。上式より融点を超える場合（\(T>T_{m}\)）、液相が安定(\(g_{S}>g_{L}\)）となり、融点を下回る場合（\(T<T_{m}\)）、固相が安定(\(g_{S}<g_{L}\)）となることが分かります。
式\((14)\)を式\((8)\)に代入することで以下の温度場を考慮したアレン-カーン方程式を得ることが出来ます。
\frac{\partial \phi}{\partial t}&=-M_{\phi}\left[30l\left(\frac{T-T_{m}}{T_{m}}\right)\phi^2\left(1-\phi\right)^2+2W\phi (1-\phi)(1-2\phi)
\right.\\
&\left. \quad\quad\quad\quad\quad-2a\left(\frac{\partial a}{\partial x} \frac{\partial \phi }{\partial x}+\frac{\partial a}{\partial y} \frac{\partial \phi }{\partial y} \right)-a^2\left(\frac{\partial ^2\phi}{\partial x^2}+\frac{\partial ^2\phi}{\partial y^2} \right)
\right.\\
&\left.\quad\quad\quad\quad\quad\quad
-\frac{\partial }{\partial x}\left( a\bar{a}k\xi \sin \left\{ k\left( \theta -\theta _0\right)\right\} \frac{\partial \phi }{\partial y}
\right)+
\frac{\partial }{\partial y}
\left(
a\bar{a}k\xi \sin \left\{ k\left( \theta -\theta _0\right)\right\} \frac{\partial \phi }{\partial x}
\right)
\right]
\tag{15}
\end{align}
上式より、フェーズフィールド変数の時間発展が温度の影響を受けるようになります。
熱伝導方程式へのフェーズフィールド変数の導入
前節では、温度場を考慮したフェーズフィールド変数の時間発展式\((15)\)を得ました。一方でフェーズフィールド変数を考慮した温度場の時間発展式も必要となります。一般に温度場の時間発展は以下の非定常熱伝導方程式によって解かれます。
c \rho\frac{\partial T}{\partial t} = K\left(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}\right) + Q\tag{16}
\end{align}
ここで\(c\)は比熱\(\rm{[J/kg\cdot K]}\)、\(\rho\)は密度\(\rm{[kg/m^3]}\)、\(K\)は熱伝導率\(\rm{[W/m \cdot K]}\)、\(Q\)は内部発熱量\(\rm{[W/m^3]}\)を表します。
また、単位体積当たりの熱容量\( c_{v}=c \rho\rm{[J/K\cdot m^3]}\)を用いて、次のように表記される場合もあります。
c_{v}\frac{\partial T}{\partial t} = K\left(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}\right) + Q\tag{17}
\end{align}
フェーズフィールド法では、内部発熱量\(Q\)の中で次のように潜熱を表現します。
Q &= l\frac{\partial h(\phi)}{\partial t}=l\frac{\partial h(\phi)}{\partial \phi}\frac{\partial \phi}{\partial t}\tag{18}\\
h(\phi) &=p(\phi)= \phi ^3(10-15\phi+6\phi^2)\tag{19}\\
\frac{\partial h(\phi)}{\partial \phi} &=30\phi^2(1-\phi)^2\tag{20}
\end{align}
ここで\(h(\phi)\)は化学的自由エネルギー密度\(g_{chem}= p(\phi)g_{a}+(1-p(\phi))g_{b}\)の定義に用いた\( p(\phi)\)と同じ関数です。グラフ化すると以下のようになります。
式\((18)\)の潜熱項は最初はイメージが付きづらいかもしれませんが、実際はそれほど難しくはありません。具体的に考えてみましょう。
例えば、\( t = 0\rightarrow 1\)の間に\(\phi=0\rightarrow 0.3\)に変化したとします。これは完全な液相が少しだけ固相になったことに相当します。この時、\(\frac{\partial \phi}{\partial t}=\frac{0.3-0}{1-0}=0.3\)です。したがって、単純には時間当たりの発熱量\(Q\)は\(0.3l\rm{[W/m^3]}\)となるはずです。（相変化が3割進むのだから潜熱も3割だけ発生するはず。）
一方で、式\((18)\)より実際の\(Q\)は次のように計算されます。
Q &=l\frac{\partial h(\phi)}{\partial \phi}\frac{\partial \phi}{\partial t}\\
&=l \times \frac{h(0.3)-h(0)}{0.3-0} \times 0.3\\
&=l \times \frac{0.163}{0.3} \times 0.3\\
&=0.163l\tag{21}
\end{align}
つまり実際の\(Q\)は\(3l\)より小さくなります。これは\(\phi\)が小さい範囲では、\(h(\phi)\)が単純な線形より緩やかに変化するためです。
すなわち、\(h(\phi)\)を導入することで、相変化の最初と最後は緩やかな発熱、途中は急速に発熱するような潜熱モデルを表現することが出来るわけです。相変化に対し単純にリニアに潜熱が発生するよりも、こちらのほうがより実際の現象に近いということなのでしょう。
以上、式\((17)\sim(20)\)をまとめると、相変化による潜熱を考慮した熱伝導方程式は次式で与えられます。
c_{v}\frac{\partial T}{\partial t} = K\left(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}\right) + 30l\phi^2 (1-\phi)^2\frac{\partial \phi}{\partial t}\tag{22}
\end{align}
これにより、熱伝導方程式にフェーズフィールド変数が導入されました。実際には式\((15)\)と式\((22)\)を連成し、交互に解くことで時間発展させます。
物性値との関連付け
実際の計算では、計算パラメーターは物性値で与えられることが多いです。そこで方程式上の係数(勾配エネルギー係数\(\bar{a}\)、ダブルウェルポテンシャル\(W\)、フェーズフィールドモビリティ\(M_{\phi}\)）と物性値(界面エネルギー\(\gamma\)、界面モビリティ\(M\)、界面幅\(\delta\))の関連付けをしておきます。
天下りではありますが、関連付けの結果は以下のようになります。
\bar{a} &= \sqrt{\frac{3\delta \gamma}{b}}
\tag{23}\\
W&=\frac{6\gamma b}{\delta}
\tag{24}\\
M_{\phi}&=\frac{bT_{m}\mu}{3\delta L}
\tag{25}\\
\end{align}
ここで\(b\)は界面幅\(\delta\)を定義する際に現れるパラメーターで\(b=2\tanh^{-1}\left(1-2\lambda\right)\)です。例えば\(\lambda=0.1\)とすると\(b \sim 2.2\)になります。
\(\bar{a},W\)（式\((23),(24)\)）に関しては過去の記事で求めたもの(界面異方性と温度場を未考慮の場合)と同等の形をしていますので導出は省略します。
\(M_{\phi}\)（式\((25)\)）に関しては過去の記事とは少々異なり、界面カイネティック係数\(\mu\rm{[m/(K \cdot s)]}\)と呼ばれるパラメーターで表現されています。この導出は特に難しいものではないのですが、以下のトグルに記載しておきますので、興味あればご覧ください。
二つの相a,bの化学的自由エネルギー\(g_{a},g_{b}\)に差があるとき、界面の移動速度\(V \rm{[m/s]}\)は、界面モビリティ\(M \rm{[m^3/J\cdot s]}\)によって次のように特徴づけられました。（過去の記事参照）
V=M\Delta g
\tag{A1}
\end{align}
ここで\(\Delta g = g_{a}-g_{b}\)は化学的駆動力（2つの化学的自由エネルギー密度の差）を表します。
この表現は粒成長や相変態問題においてよく用いられます。
一方で凝固問題の場合、下式のような融点を用いた表現が用いられることが多いようです。
V=\mu \left( T_m-T\right)
\tag{A2}
\end{align}
ここで\(\mu\)は界面カイネティック係数\(\rm{[m/(K \cdot s)]}\)と呼ばれる正の定数です。
上式は、界面エネルギーを融点で表現しており、物質の温度\(T\)が融点\(T_{m}\)を上回るか・下回るかで界面の移動方向が決まります。例えば\(T>T_{m}\)の時、固相が液相になる方向へ界面が移動します。
さて、界面速度\(V\)はアレン-カーン方程式の変形からも得ることが出来ました。ベーシックな平衡プロファイルを有する界面が、一定速度\(V\)で\(+x\)方向に移動する定常成長問題では次式が成り立ちます。（過去の記事:式\((37)\)参照）
V=\frac{6a}{\sqrt{2W}}M_{\phi}\Delta g \tag{A3}
\end{align}
上式に式\((8)\)を代入すると、
V&=-\frac{6a}{\sqrt{2W}}M_{\phi}l\left(\frac{T-T_{m}}{T_{m}}\right)\\
&=\frac{6aM_{\phi}l}{\sqrt{2W}T_{m}}\left(T_{m}-T\right) \tag{A4}
\end{align}
上式を式\((A2)\)と比較すると、次式が成り立ちます。
&\mu \left( T_m-T\right)=\frac{6aM_{\phi}l}{\sqrt{2W}T_{m}}\left(T_{m}-T\right) \\
&\rightarrow\mu =\frac{6aM_{\phi}l}{\sqrt{2W}T_{m}} \\
&\rightarrow M_{\phi}=\frac{T_{m}\mu} {6al}\sqrt{2W}
\tag{A5}
\end{align}
式\((17),(18)\)を代入すると
M_{\phi}&=\frac{T_{m}\mu} {6l\sqrt{\frac{3\delta \gamma}{b}}}\sqrt{2\cdot\frac{6\gamma b}{\delta}}\\
&=\frac{T_{m}\mu} {6l}\sqrt{\frac{b}{3\delta \gamma}\cdot\frac{12\gamma b}{\delta}}\\
&=\frac{bT_{m}\mu}{3\delta L}
\tag{A6}
\end{align}
以上で式\((15)\)が導出できました。
おわりに
今回は純物質凝固のモデル化についてまとめました。純物質という仮定により、液相と固相の化学的自由エネルギーの差が潜熱分のみとなり、駆動力がシンプルに書けるのがポイントです。次回は今回得られた方程式を離散化し、実際にシミュレーションしてみたいと思います。
参考文献
[1] 高木 知弘ら, “フェーズフィールド法”, 養賢堂, 2012/3/2
[2] 山中 晃徳ら, “Pythonによるフェーズフィールド法入門 基礎理論からデータ同化の実装まで”, 丸善出版, 2023/12/15
コメント