フェーズフィールド法⑮：WBMモデル②

はじめに

　前回はWBMモデルにおけるアレン-カーン方程式を導出を行いました。今回は濃度場の時間発展を記述するカーン-ヒリアード方程式の導出を行います。

はじめに　２元系合金の凝固に関する基本的なモデルとしては、1991年頃にWheeler、Boettinger、McFaddenらによって提案されたWBMモデルが知られています。WBMモデルはアレン-カーン方程式とカーン-ヒリアード方程式を...

濃度場に関するカーン-ヒリアード方程式の導出

　WBMモデルにおいて、濃度場の時間発展は以下のエントロピー\(S\)に関するカーン-ヒリアード方程式によって記述されます。

\begin{align}
\frac{dc}{dt}=-\nabla \cdot \left(M_{c}\nabla\frac{\delta S}{\delta c}\right)\tag{1}
\end{align}

ここで\(c\)は元素Ｂのモル分率です。また、右辺の括弧内はベクトルとなります。

上式はエントロピーが増大する方向に濃度場が進展することを示しています。

フェーズフィールド変数は非保存量であるため、アレン-カーン方程式を用いましたが、濃度は保存量であるためカーン-ヒリアード方程式を用いることに注意しましょう。

この時エントロピー\(S\)はバルクのエントロピーと界面のエントロピーの和として表され、次のように与えられました。（前回記事参照）

\begin{align}
S=\int s_{dens}dV =\int s(\phi ,e , c)-\frac{\epsilon}{2} \left| \nabla \phi\right|^2dV \tag{2}
\end{align}

従って、式\((2)\)を用いると式\((1)\)右辺の汎関数微分は以下のように書けます。

\begin{align}
\frac{\delta S}{\delta c}&=\frac{\partial s_{dens}}{\partial c}
\\
&=\frac{\partial}{\partial c}\left( s(\phi ,e , c)-\frac{\epsilon ^2}{2} \left| \nabla \phi\right|^2\right)\\
&=\frac{\partial s(\phi ,e , c)}{\partial c} \tag{3}
\end{align}

また、\(\frac{\partial s(\phi ,e , c)}{\partial c} \)は化学ポテンシャル\(\mu^A,\mu^B\)を用いて、次のように表すことが出来ます。

\begin{align}
\frac{\partial s(\phi ,e , c)}{\partial c}=\frac{1}{T}\left[\left\{f_A(T,\phi)+\frac{RT}{v_{m}}\log (1-c)\right\}-\left\{f_B(T,\phi)+\frac{RT}{v_{m}}\log c\right\}\right]
\tag{4}
\end{align}

式\((4)\)の導出

　各元素の化学ポテンシャル\(\mu^A,\mu^B\)は1モル当たりのヘルムホルツの自由エネルギーと考えられるため、全自由エネルギー\(f\)は次のように書けます。

\begin{align}
f&=(1-c)\mu^A+c\mu^B\tag{A1}
\end{align}

また、正則容体近似により化学ポテンシャルは次式にて与えられます。

\begin{align}
\mu^A=f_A(T,\phi)+\Omega(\phi) c^2+\frac{RT}{v_{m}}\log (1-c)\tag{A2}\\
\mu^B=f_B(T,\phi)+\Omega(\phi) (1-c)^2+\frac{RT}{v_{m}}\log c\tag{A3}
\end{align}

したがって上式を式\((A1)\)に代入すると

\begin{align}
f&=(1-c)\left(f_A(T,\phi)+\Omega(\phi) c^2+\frac{RT}{v_{m}}\log (1-c)\right)\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+c\left(f_B(T,\phi)+\Omega(\phi) (1-c)^2+\frac{RT}{v_{m}}\log c\right)\tag{A4}
\end{align}

これを用いると、一般的な熱力学上の関係式から、エントロピー密度\(s\)は次式で表されます。（前回記事参照）

\begin{align}
s&=\frac{1}{T}(e-f)\\
&=\frac{1}{T}\left[e-(1-c)\left(f_A(T,\phi)+\Omega(\phi) c^2+\frac{RT}{v_{m}}\log (1-c)\right)\right.\\
&\left.\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad-c\left(f_B(T,\phi)+\Omega(\phi) (1-c)^2+\frac{RT}{v_{m}}\log c\right)\right]\tag{A5}
\end{align}

さらに上式を濃度\(c\)で微分すると、

\begin{align}
\frac{\partial s}{\partial c}&=\frac{1}{T}\frac{\partial}{\partial c}\left[e-(1-c)\left\{f_A(T,\phi)+\Omega(\phi) c^2+\frac{RT}{v_{m}}\log (1-c)\right\}\right.\\
&\left.\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad-c\left\{f_B(T,\phi)+\Omega(\phi) (1-c)^2+\frac{RT}{v_{m}}\log c\right\}\right]\tag{A6}
\end{align}

ここでWBMモデルでは、温度\(T\)および内部エネルギ\(e\)は一定、かつ混合エンタルピー\(\Omega=0\)を仮定するので（前回記事参照）、上式は次のように変形できます。

\begin{align}
\frac{\partial s}{\partial c}&=\frac{1}{T}\frac{\partial}{\partial c}\left[-(1-c)\left\{f_A(T,\phi)+\frac{RT}{v_{m}}\log (1-c)\right\}-c\left\{f_B(T,\phi)+\frac{RT}{v_{m}}\log c\right\}\right]\\
&=\frac{1}{T}\left[\left\{f_A(T,\phi)+\frac{RT}{v_{m}}\log (1-c)\right\}-(1-c)\left\{\frac{RT}{v_{m}}\frac{-1}{1-c}\right\}\right.\\
&\left.\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad-\left\{f_B(T,\phi)+\frac{RT}{v_{m}}\log c\right\}-c\left\{\frac{RT}{v_{m}}\frac{1}{c}\right\}\right]\\
&=\frac{1}{T}\left[\left\{f_A(T,\phi)+\frac{RT}{v_{m}}\log (1-c)\right\}+\frac{RT}{v_{m}}-\left\{f_B(T,\phi)+\frac{RT}{v_{m}}\log c\right\}-\frac{RT}{v_{m}}\right]\\
&=\frac{1}{T}\left[\left\{f_A(T,\phi)+\frac{RT}{v_{m}}\log (1-c)\right\}-\left\{f_B(T,\phi)+\frac{RT}{v_{m}}\log c\right\}\right]
\tag{A7}
\end{align}

以上で式\((4)\)が求まりました。

ちなみに式\((A7)\)は\(\Omega=0\)を踏まえた正則容体近似式\((A2),(A3)\)を用いると以下のように書くこともできます。

\begin{align}
\frac{\partial s}{\partial c}=\frac
{1}{T}\left(\mu^A-\mu^B\right)\tag{A8}
\end{align}

よってカーン-ヒリアード方程式\((1)\)は次のように書けます。

\begin{align}
\frac{dc}{dt} &= -\nabla \cdot \left[M_{c}\nabla \frac{1}{T} \left\{ \left(f_A + \frac{RT}{v_{m}} \log (1-c) \right) – \left(f_B + \frac{RT}{v_{m}} \log c \right) \right\} \right] \tag{5}
\end{align}

上式右辺は以下のように変形できます。

\begin{align}
-\nabla \cdot &\left[M_{c}\nabla \frac{1}{T} \left\{ \left(f_A + \frac{RT}{v_{m}} \log (1-c) \right) – \left(f_B + \frac{RT}{v_{m}} \log c \right) \right\} \right]\\
&=-\nabla \cdot \left[M_{c}\nabla\left\{\left(\frac{f_A}{T}+\frac{R}{v_{m}}\log (1-c)\right)-\left(\frac{f_B}{T}+\frac{R}{v_{m}}\log c\right) \right\}\right]\\
&=-\nabla \cdot \left[M_{c}\left\{\left(\frac{1}{T}\nabla f_A+\frac{R}{v_{m}}\nabla\log (1-c)\right)-\left(\frac{1}{T}\nabla f_B+\frac{R}{v_{m}}\nabla\log c\right) \right\}\right]\\

&=-\nabla \cdot \left[M_{c}\left\{\left(\frac{1}{T}
\begin{pmatrix}
\frac{\partial f_A}{\partial x}\\
\frac{\partial f_A}{\partial y}\\
\end{pmatrix}
+\frac{R}{v_{m}}
\begin{pmatrix}
\frac{\partial \log (1-c)}{\partial x}\\
\frac{\partial \log (1-c)}{\partial y}\\
\end{pmatrix}
\right)-\left(\frac{1}{T}
\begin{pmatrix}
\frac{\partial f_B}{\partial x}\\
\frac{\partial f_B}{\partial y}\\
\end{pmatrix}
+\frac{R}{v_{m}}
\begin{pmatrix}
\frac{\partial \log c}{\partial x}\\
\frac{\partial \log c}{\partial y}\\
\end{pmatrix}
\right) \right\}\right]\\
&=-\nabla \cdot \left[M_{c}\left\{\left(\frac{1}{T}
\begin{pmatrix}
\frac{\partial f_A}{\partial \phi}\frac{\partial \phi}{\partial x}\\
\frac{\partial f_A}{\partial \phi}\frac{\partial \phi}{\partial y}\\
\end{pmatrix}
+\frac{R}{v_{m}}
\begin{pmatrix}
\frac{\partial \log (1-c)}{\partial (1-c)}\frac{\partial \log (1-c)}{\partial c}\frac{\partial c}{\partial x}\\
\frac{\partial \log (1-c)}{\partial (1-c)}\frac{\partial \log (1-c)}{\partial c}\frac{\partial c}{\partial y}\\
\end{pmatrix}\right)\right.\right.\\
&\left.\left.\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad
-\left(\frac{1}{T}
\begin{pmatrix}
\frac{\partial f_B}{\partial \phi}\frac{\partial \phi}{\partial x}\\
\frac{\partial f_B}{\partial \phi}\frac{\partial \phi}{\partial y}\\
\end{pmatrix}
+\frac{R}{v_{m}}
\begin{pmatrix}
\frac{\partial \log c}{\partial c}\frac{\partial c}{\partial x}\\
\frac{\partial \log c}{\partial c}\frac{\partial c}{\partial y}\\
\end{pmatrix}
\right) \right\}\right]\\
&=-\nabla \cdot \left[M_{c}\left\{\left(\frac{1}{T}\frac{\partial f_A}{\partial \phi}\nabla \phi
-\frac{R}{v_{m}}\frac{1}{1-c}\nabla c\right)-\left(\frac{1}{T}\frac{\partial f_B}{\partial \phi}\nabla \phi+\frac{R}{v_{m}}\frac{1}{c}\nabla c\right) \right\}\right]\\
&=-\nabla \cdot \left[M_{c}\left\{\left(\frac{1}{T}\frac{\partial f_A}{\partial \phi}-\frac{1}{T}\frac{\partial f_B}{\partial \phi}\right)\nabla \phi-\frac{R}{v_{m}}\left(\frac{1}{1-c}+\frac{1}{c}\right) \nabla c\right\}\right]\\
&=\nabla \cdot \left[M_{c}\left\{\left(\frac{1}{T}\frac{\partial f_B}{\partial \phi}-\frac{1}{T}\frac{\partial f_A}{\partial \phi}\right)\nabla \phi+\frac{R}{v_{m}}\left(\frac{1}{c(1-c)}\right) \nabla c\right\}\right]\\
&=\nabla \cdot \left[M_{c}\left\{\left(H^B-H^A\right)\nabla \phi+\frac{R}{v_{m}}\frac{\nabla c}{c(1-c)} \right\}\right]\\
&=\nabla \cdot \left[\frac{R}{v_{m}}\frac{1}{c(1-c)}M_{c} \left\{\nabla c+\frac{v_{m}}{R}c(1-c)\left(H^B-H^A\right)\nabla \phi\right\}\right]\\
&=\nabla \cdot \left[D\left\{\nabla c+\frac{v_{m}}{R}c(1-c)\left(H^B-H^A\right)\nabla \phi\right\}\right]
\tag{6}
\end{align}

ここで\(H^{A},H^{B}\)は前回記事同様、下記を表しています。

\begin{align}
H^{A}&=\frac{1}{T}\frac{\partial f_A}{\partial \phi}=W^{A}\frac{\partial g}{\partial \phi}+30gL^{A} \left(\frac{1}{T}-\frac{1}{T_{m}^{A}}\right)\tag{7}\\
H^{B}&=\frac{1}{T}\frac{\partial f_B}{\partial \phi}=W^{B}\frac{\partial g}{\partial \phi}+30gL^{B} \left(\frac{1}{T}-\frac{1}{T_{m}^{B}}\right)\tag{8}
\end{align}

また\(D\)は次のように置きました。

\begin{align}
D&=\frac{R}{v_{m}}\frac{1}{c(1-c)}M_{c} \tag{9}
\end{align}

したがって、式\((6)\)を式\((5)\)右辺に代入すると、濃度場に関するカーン-ヒリアード方程式は最終的に以下のように表せます。

\begin{align}
\frac{dc}{dt}&=\nabla \cdot \left[D\left\{\nabla c+\frac{v_{m}}{R}c(1-c)\left(H^B-H^A\right)\nabla \phi\right\}\right]
\tag{10}
\end{align}

拡散係数

　式\((10)\)がどのような性質の解をもつのかパッと見分かりづらいかもしれません。

例えば、固相（\(\phi = 0\)）または液相（\(\phi = 1\)）領域を仮定すると、\(\nabla \phi= 0\)となり、式\((10)\)は次のように変形できます。

\begin{align}
\frac{dc}{dt}=\nabla \cdot D \nabla c\tag{11}
\end{align}

仮に\(D\)が定数と考えると、上式はいわゆる拡散方程式そのものです。つまり\(D\)は拡散係数に相当するパラメーターとなります。

このことは仮に固相・液相中に濃度ムラがあった場合、それが均一になるように濃度場が進展することを示しています。

また今回は凝固現象を取り扱うため、\(D\)は固体・液体の拡散係数\(D_{S},D_{L}\)を用いて次のように書くこととします。

\begin{align}
D=D_{S}+p(\phi)(D_L-D_S)\tag{12}
\end{align}

つまり相分率\(\phi\)を重みとした固相と液相の重ね合わせで、系全体の拡散係数を決めるわけです。

このように書くと単一相中の拡散係数は\(D_S\)または\(D_L\)のいずれかの定数となり、式\((10)\)は拡散方程式に帰着します。一方で固液界面（\(\phi \neq 0,1\)）では相分率に応じた拡散係数となります。

上式と式\((9)\)と比較すると、\(M_c\)は次のように書くことが出来ます。

\begin{align}
D_{S}&+p(\phi)(D_L-D_S)=\frac{R}{v_{m}}\frac{1}{c(1-c)}M_{c} \\
&\rightarrow M_{c}=\frac{v_{m}}{R}c(1-c)\left[D_{S}+p(\phi)(D_L-D_S)\right] \tag{13}
\end{align}

一般に\(D_S \sim 10^{-13}\rm{[m^2/s]}\)、\(D_L \sim 10^{-9}\rm{[m^2/s]}\)程度となり、固相より液相の拡散係数のほうが、\(10^4\)倍程大きくなることが知られています。

おわりに

　今回はWBMモデルのカーン-ヒリアード方程式の導出を行いました。以上で一通り支配方程式が出そろいましたので、次回は1次元解析で方程式の解の性質を調べていきたいと思います。

参考文献

[1] 高木知弘ら, “フェーズフィールド法”, 養賢堂, 2012/3/2
[2] 山中晃徳ら, “Pythonによるフェーズフィールド法入門基礎理論からデータ同化の実装まで”, 丸善出版, 2023/12/15
[3] 小山敏幸ら, “フェーズフィールド法入門 (計算力学レクチャーコース) “, 丸善出版, 2013/4/13
[4] J.A. Warren, W.J. Boettinger : “Prediction of dendritic growth and microsegregation patterns in a binary alloy using the phase-field method”. Acta Metallurgica et Materialia, Volume 43, Issue 2, February 1995, Pages 689-703
[5] S.-L. Wang, R.F. Sekerka : “Thermodynamically-consistent phase-field models for solidification”. Physica D: Nonlinear Phenomena, Volume 69, Issues 1-2, 15 November 1993, Pages 189-200