はじめに
前回はデンドライト構造形成のメカニズムについてまとめました。デンドライト構造は過冷却状態において、化学的駆動力のゆらぎによって生じた微小な突起構造が熱的に安定化しようと発達した結果生じる結晶構造でした。そこで今回は化学的駆動力のゆらぎを定式化し、純物質におけるデンドライト形成のシミュレーション行います。
化学的駆動力のゆらぎのモデル化
マクロスケールでは温度が一定に見えたとしても、ミクロスケールでは分子や原子の運動によるエネルギーのゆらぎが存在します。これが局所的な化学的自由エネルギーの変動、すなわち化学的駆動力のゆらぎを引き起こします。
化学的駆動力のゆらぎ\(\varepsilon\)は次のように定式化されます[1]。
\begin{align}
\varepsilon=4W\phi\left(1-\phi\right)\chi\tag{1}\\
\end{align}
\varepsilon=4W\phi\left(1-\phi\right)\chi\tag{1}\\
\end{align}
ここで\(W\)はダブルウェルポテンシャルの大きさを決める係数、\(\chi\)は乱数で\(-0.1\leq \chi \leq 0.1\)の間の値を取るものとします。また、係数の4は文献[1]において後の式変形をスマートにするためにかけられており、物理的な意味はありません。
このままだと上式の意味が掴みづらいので、\(\chi\)の係数部分\(4W\phi\left(1-\phi\right)\)について、いくつかの\(W\)でグラフ化してみます。
グラフより\(\phi=0.5\)、すなわち固液界面で最も値が大きくなることが分かります。これは界面のゆらぎが最も大きく、単一相に近づくにつれて小さくなっていくことを示しています。
また、\(W\)が大きくなるにつれて係数の値も大きくなっています。\(W\)はダブルウェルポテンシャルにおける障壁の高さを表していました。すなわち異なる相に遷移しようとするとき、ギブスのエネルギー的に高い状態となり、それに応じて界面の不安定性も増加することを示しています。
以上より化学的駆動力のゆらぎ\(\varepsilon\)の影響は、単純な乱数で与えられるわけではなく、固液界面に近づくほど大きく、また、その相変化にかかるエネルギー（障壁高さ）が大きいほど大きくなることが分かりました。
ゆらぎを考慮したアレン-カーン方程式
次にゆらぎの効果（式\((1)\)）を実際にアレン-カーン方程式に組み込みます。アレン-カーン方程式はもともとギブスの自由エネルギーの汎関数微分として、次のように与えられました。（過去記事参照）
\begin{align}
\frac{\partial \phi}{\partial t}&=-M_{\phi}\frac{\delta G}{\delta \phi} \tag{2}\\
\end{align}
\frac{\partial \phi}{\partial t}&=-M_{\phi}\frac{\delta G}{\delta \phi} \tag{2}\\
\end{align}
上式に対して、ゆらぎの効果（式\((1)\)）は次式のように加えられます。
\begin{align}
\frac{\partial \phi}{\partial t}&=-M_{\phi}\left(\frac{\delta G}{\delta \phi}+\varepsilon\right) \tag{3}\\
\end{align}
\frac{\partial \phi}{\partial t}&=-M_{\phi}\left(\frac{\delta G}{\delta \phi}+\varepsilon\right) \tag{3}\\
\end{align}
自由エネルギーにランダムなゆらぎを加えればいいので、必ずしも上式のような形である必要はない気もしますが、ここでは文献[1]に倣った定式化をすることとします。
さらに上式を純物質凝固のモデルとして変形すると、
\begin{align}
\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-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
-\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)
+4W\phi\left(1-\phi\right)\chi\right] \tag{4}\\
\end{align}
\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-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
-\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)
+4W\phi\left(1-\phi\right)\chi\right] \tag{4}\\
\end{align}
上式のままでもいいのですが、文献[1]に倣った式変形をしていくと、以下のように続きます。
\begin{align}
\frac{\partial \phi}{\partial t}&=-M_{\phi}\left[4W\phi\left(1-\phi\right)\left\{\frac{15}{2W}l\left(\frac{T-T_{m}}{T_{m}}\right)\phi\left(1-\phi\right)+\frac{1}{2}(1-2\phi)+\chi\right\}
\right.\\
&\left. \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
-\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] \\
\rightarrow
&\frac{\partial \phi}{\partial t}=M_{\phi}\left[4W\phi\left(\phi-1\right)\left\{\frac{15}{2W}l\left(\frac{T-T_{m}}{T_{m}}\right)\phi\left(1-\phi\right)+\phi-0.5+\chi\right\}
\right.\\
&\left. \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
+\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{5}\\
\end{align}
\frac{\partial \phi}{\partial t}&=-M_{\phi}\left[4W\phi\left(1-\phi\right)\left\{\frac{15}{2W}l\left(\frac{T-T_{m}}{T_{m}}\right)\phi\left(1-\phi\right)+\frac{1}{2}(1-2\phi)+\chi\right\}
\right.\\
&\left. \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
-\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] \\
\rightarrow
&\frac{\partial \phi}{\partial t}=M_{\phi}\left[4W\phi\left(\phi-1\right)\left\{\frac{15}{2W}l\left(\frac{T-T_{m}}{T_{m}}\right)\phi\left(1-\phi\right)+\phi-0.5+\chi\right\}
\right.\\
&\left. \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
+\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{5}\\
\end{align}
最後の式変形には\(\chi=-\chi\)を利用しています。以上がゆらぎを導入した純物質凝固のアレン-カーン方程式です。
さらに式\((4)\)について離散化します。純物凝固モデルにゆらぎを足せばいいだけですので簡単です。
\begin{align}
\frac{\phi_{i,j}^{n+1}-\phi_{i,j}^{n}}{\Delta t}&=-M_{\phi}\left[6l\left(\frac{T_{i,j}^{n}-T_{m}}{T_{m}}\right)\phi_{i,j}^{n}\left(1-\phi_{i,j}^{n}\right)+2W\phi_{i,j}^{n} (1-\phi_{i,j}^{n})(1-2\phi_{i,j}^{n})
\right.\\
&\left. -2a_{i,j}^{n}\left(\frac{a_{i+1,j}^{n}-a_{i-1,j}^{n}}{2\Delta x} \frac{ \phi_{i+1,j}^{n}-\phi_{i-1,j}^{n} }{2\Delta x}+\frac{a_{i,j+1}^{n}-a_{i,j-1}^{n}}{2\Delta y} \frac{\phi_{i,j+1}^{n}-\phi_{i,j-1}^{n} }{2\Delta y} \right)\\
\quad\quad-\left(a_{i,j}^{n}\right)^2\Delta^9\phi_{i,j}^{n}
-\frac{A_{i+1,j}^{n}-A_{i-1,j}^{n} }{2\Delta x}+
\frac{B_{i,j+1}^{n}-B_{i,j-1}^{n} }{2\Delta y} +4W\phi_{i,j}\left(1-\phi_{i,j}\right)\chi
\right] \\
\rightarrow
\phi_{i,j}^{n+1}&=\phi_{i,j}^{n}+M_{\phi}\Delta t\left[-6l\left(\frac{T_{i,j}^{n}-T_{m}}{T_{m}}\right)\phi_{i,j}^{n}\left(1-\phi_{i,j}^{n}\right)-2W\phi_{i,j}^{n} (1-\phi_{i,j}^{n})(1-2\phi_{i,j}^{n})
\right.\\
&\left. +2a_{i,j}^{n}\left(\frac{a_{i+1,j}^{n}-a_{i-1,j}^{n}}{2\Delta x} \frac{ \phi_{i+1,j}^{n}-\phi_{i-1,j}^{n} }{2\Delta x}+\frac{a_{i,j+1}^{n}-a_{i,j-1}^{n}}{2\Delta y} \frac{\phi_{i,j+1}^{n}-\phi_{i,j-1}^{n} }{2\Delta y} \right)\\
\quad\quad+\left(a_{i,j}^{n}\right)^2\Delta^9\phi_{i,j}^{n}+\frac{A_{i+1,j}^{n}-A_{i-1,j}^{n} }{2\Delta x}-
\frac{B_{i,j+1}^{n}-B_{i,j-1}^{n} }{2\Delta y} +4W\phi_{i,j}\left(1-\phi_{i,j}\right)\chi
\right] \tag{6}\\
A_{i,j}^{n}&=a_{i,j}^{n}\bar{a}k\xi \sin \left\{ k\left( \theta_{i,j}^{n} -\theta _0\right)\right\} \frac{\phi_{i,j+1}^{n}-\phi_{i,j-1}^{n} }{2\Delta y} \tag{7}\\
B_{i,j}^{n}&=a_{i,j}^{n}\bar{a}k\xi \sin \left\{ k\left( \theta_{i,j}^{n} -\theta _0\right)\right\} \frac{\phi_{i+1,j}^{n}-\phi_{i-1,j}^{n} }{2\Delta x}\tag{8}\\
a_{i,j}^{n}&=\bar{a}\left[1+\xi \cos \left\{ k\left( \theta_{i,j}^{n} -\theta _0\right)\right\} \right]
\tag{9}
\end{align}
\frac{\phi_{i,j}^{n+1}-\phi_{i,j}^{n}}{\Delta t}&=-M_{\phi}\left[6l\left(\frac{T_{i,j}^{n}-T_{m}}{T_{m}}\right)\phi_{i,j}^{n}\left(1-\phi_{i,j}^{n}\right)+2W\phi_{i,j}^{n} (1-\phi_{i,j}^{n})(1-2\phi_{i,j}^{n})
\right.\\
&\left. -2a_{i,j}^{n}\left(\frac{a_{i+1,j}^{n}-a_{i-1,j}^{n}}{2\Delta x} \frac{ \phi_{i+1,j}^{n}-\phi_{i-1,j}^{n} }{2\Delta x}+\frac{a_{i,j+1}^{n}-a_{i,j-1}^{n}}{2\Delta y} \frac{\phi_{i,j+1}^{n}-\phi_{i,j-1}^{n} }{2\Delta y} \right)\\
\quad\quad-\left(a_{i,j}^{n}\right)^2\Delta^9\phi_{i,j}^{n}
-\frac{A_{i+1,j}^{n}-A_{i-1,j}^{n} }{2\Delta x}+
\frac{B_{i,j+1}^{n}-B_{i,j-1}^{n} }{2\Delta y} +4W\phi_{i,j}\left(1-\phi_{i,j}\right)\chi
\right] \\
\rightarrow
\phi_{i,j}^{n+1}&=\phi_{i,j}^{n}+M_{\phi}\Delta t\left[-6l\left(\frac{T_{i,j}^{n}-T_{m}}{T_{m}}\right)\phi_{i,j}^{n}\left(1-\phi_{i,j}^{n}\right)-2W\phi_{i,j}^{n} (1-\phi_{i,j}^{n})(1-2\phi_{i,j}^{n})
\right.\\
&\left. +2a_{i,j}^{n}\left(\frac{a_{i+1,j}^{n}-a_{i-1,j}^{n}}{2\Delta x} \frac{ \phi_{i+1,j}^{n}-\phi_{i-1,j}^{n} }{2\Delta x}+\frac{a_{i,j+1}^{n}-a_{i,j-1}^{n}}{2\Delta y} \frac{\phi_{i,j+1}^{n}-\phi_{i,j-1}^{n} }{2\Delta y} \right)\\
\quad\quad+\left(a_{i,j}^{n}\right)^2\Delta^9\phi_{i,j}^{n}+\frac{A_{i+1,j}^{n}-A_{i-1,j}^{n} }{2\Delta x}-
\frac{B_{i,j+1}^{n}-B_{i,j-1}^{n} }{2\Delta y} +4W\phi_{i,j}\left(1-\phi_{i,j}\right)\chi
\right] \tag{6}\\
A_{i,j}^{n}&=a_{i,j}^{n}\bar{a}k\xi \sin \left\{ k\left( \theta_{i,j}^{n} -\theta _0\right)\right\} \frac{\phi_{i,j+1}^{n}-\phi_{i,j-1}^{n} }{2\Delta y} \tag{7}\\
B_{i,j}^{n}&=a_{i,j}^{n}\bar{a}k\xi \sin \left\{ k\left( \theta_{i,j}^{n} -\theta _0\right)\right\} \frac{\phi_{i+1,j}^{n}-\phi_{i-1,j}^{n} }{2\Delta x}\tag{8}\\
a_{i,j}^{n}&=\bar{a}\left[1+\xi \cos \left\{ k\left( \theta_{i,j}^{n} -\theta _0\right)\right\} \right]
\tag{9}
\end{align}
以上が離散化形式です。あとは前回記事同様、上式と潜熱を考慮した非定常熱伝導方程式を連成して解けばOKです。
解析条件
では実際にシミュレーションしてみます。解析対象は過去の記事と同様に過冷却状態における純Niの凝固とします。Niの物性は次の通りです。
\begin{align}
\begin{array}{ccc}
\hline
\textbf{Ni Property} & \textbf{value} & \textbf{units}\\
\hline
融点 T_{m} & 1728 & \rm{[K]} \\
界面エネルギー \gamma & 0.37 & \rm{[J/m^2]} \\
熱伝導率 K & 84.01 & \rm{[W/m⋅K]} \\
体積比熱 c_{v} & 5.42 \times 10^{6} & \rm{ [J/K⋅m^3]} \\
潜熱密度 l & 2.35 \times 10^{9}& \rm{[J/mol]}\\
界面カイネティック係数 \mu & 2.0 & \rm{[m/(K⋅s)]}\\
\hline
\end{array}
\end{align}
\begin{array}{ccc}
\hline
\textbf{Ni Property} & \textbf{value} & \textbf{units}\\
\hline
融点 T_{m} & 1728 & \rm{[K]} \\
界面エネルギー \gamma & 0.37 & \rm{[J/m^2]} \\
熱伝導率 K & 84.01 & \rm{[W/m⋅K]} \\
体積比熱 c_{v} & 5.42 \times 10^{6} & \rm{ [J/K⋅m^3]} \\
潜熱密度 l & 2.35 \times 10^{9}& \rm{[J/mol]}\\
界面カイネティック係数 \mu & 2.0 & \rm{[m/(K⋅s)]}\\
\hline
\end{array}
\end{align}
またメッシュ数等、その他計算パラメーターは下表のとおり設定します。
\begin{align}
\begin{array}{ccc}
\hline
\textbf{calculation parameters} & \textbf{value} & \textbf{units}\\
\hline
メッシュ数 n_{x},n_{y} & 750 & \rm{[-]} \\
メッシュサイズ dx,dy & 2.0\times10^{-8} & \rm{[m]} \\
タイムステップ dt & 5.0\times10^{-12}& \rm{[s]} \\
界面幅 \delta & 4dx & \rm{[m]} \\
異方性強度 \xi & 0.03 & \rm{ [-]} \\
異方性モード k & 4 & \rm{[-]}\\
優先成長方向 \theta_0 & 0 & \rm{[deg]}\\
\hline
\end{array}
\end{align}
\begin{array}{ccc}
\hline
\textbf{calculation parameters} & \textbf{value} & \textbf{units}\\
\hline
メッシュ数 n_{x},n_{y} & 750 & \rm{[-]} \\
メッシュサイズ dx,dy & 2.0\times10^{-8} & \rm{[m]} \\
タイムステップ dt & 5.0\times10^{-12}& \rm{[s]} \\
界面幅 \delta & 4dx & \rm{[m]} \\
異方性強度 \xi & 0.03 & \rm{ [-]} \\
異方性モード k & 4 & \rm{[-]}\\
優先成長方向 \theta_0 & 0 & \rm{[deg]}\\
\hline
\end{array}
\end{align}
界面の曲率評価を正確に行うため、メッシュ数は多めに\(n_{x}=n_{y}=750\)とします。
初期条件は、過冷却状態の溶融Ni液体（\(\phi=0\)）の左下に半径\(2dx\)、温度がちょうど融点（T=1728[K]）のNi固体（\(\phi=1\)）の凝固核を配置した状態とします。
また過冷度\(\Delta =\frac{c_v\left(T_{m}-T_{0}\right)}{l} \)をパラメーターとし、溶融Ni液体の温度として\(T_{0}=1597.9/1511.2/1424.5\rm{[K]}(\Delta=0.3/0.5/0.7)\)の三条件を比較することとします。
MATLABコード
以下に今回使用したMATLABコードを記載します。
% 初期化処理 --------------------------------------------------------------
clc; close all;
clear all;
% パラメーター設定 --------------------------------------------------------
% 時間・空間ステップ設定
dt = 5.0e-12; % タイムステップ[s]
maxstep = 45000; % 最大ステップ数
dx = 2.0e-08; % 差分格子点の間隔[m]:20[um]
dy = dx; % 差分格子点の間隔[m]
nx = 750; % 差分格子点数
ny = nx; % 差分格子点数
xmax = dx * nx ;
xmin = 0;
ymax = dy * ny;
ymin = 0;
xl = xmin:dx:xmax-dx;
yl = ymin:dy:ymax-dy;
% 物性値の設定
T_0 = 1597.9; % 系の温度 [K]
T_melt = 1728.0; % Niの融点 [K]
gamma = 0.37; % 界面エネルギー [J/m2]
K = 84.01; % 熱伝導率 [W/(mK)]
cv = 5.42e+06; % 比熱 [J/K]
latent = 2.35e+09; % 潜熱 [J/mol]
mu = 2.0; % 界面カイネティック係数 [m/(Ks)]
lamb = 0.1; % 界面厚みの閾値
b = 2.0 * atanh(1.0 - 2.0 * lamb);
delta = 4.0 * dx; % 界面幅[m]
sc = cv *(T_melt-T_0)/latent; % 過冷度
% 物性値の変換
w = 6.0 * gamma * b / delta; % エネルギー障壁[J/m3]
mp = b * T_melt * mu / (3.0 * delta * latent); % フェーズフィールドモビリティー
abar = sqrt(3.0*delta*gamma/b); % 基準となる勾配係数
xi = 0.05; % 異方性強度
k = 4; % 異方性モード
theta_0 = 0; % 優先成長方向
% 変数宣言
phi = zeros(nx,ny);
phi_new = zeros(nx,ny);
dphi = zeros(nx,ny);
q = zeros(nx,ny);
nvec_x = zeros(nx,ny);
nvec_y = zeros(nx,ny);
T = zeros(nx,ny);
dT = zeros(nx,ny);
T_new = zeros(nx,ny);
dg = zeros(nx,ny);
theta = zeros(nx,ny);
thetadeg = zeros(nx,ny);
a = zeros(nx,ny);
A = zeros(nx,ny);
B = zeros(nx,ny);
dadx= zeros(nx,ny);
dady = zeros(nx,ny);
dphidx = zeros(nx,ny);
dphidy = zeros(nx,ny);
lap_phi = zeros(nx,ny);
lap_temp = zeros(nx,ny);
% 安定条件の確認
theta_test = 0:2*pi/100:2*pi;
amax = max(abar * (1 + xi * cos(k*(theta_test-theta_0))));% 勾配係数の最大値を計算する。
dt1 = (1/(2*amax^2*mp))*(dx^2*dy^2)/(dx^2+dy^2);
dt2 =(cv/(2*K))*(dx^2*dy^2)/(dx^2+dy^2);
if dt > min(dt1,dt2)
error('タイムステップが安定条件を満たしていません。')
end
% 可視化設定 -----------------------------------------------------------------
% ファイル名
paramstr =['abar = ', num2str(abar, '%.1e'), ', w = ', num2str(w, '%.1e'), ...
', dg = ', num2str(1, '%.2f'), ', mp = ', num2str(mp, '%.2f'),...
', xi = ', num2str(xi, '%.2f'), ', k = ', num2str(k, '%.2f'),...
', theta_0 = ', num2str(theta_0, '%.2f'),', Delta = ', num2str(sc, '%.2f'), ];
filename = [paramstr,'.mp4'];
% VideoWriter オブジェクトの初期化
outputVideo = VideoWriter(filename, 'MPEG-4');
outputVideo.FrameRate = 30; % フレームレートの設定
open(outputVideo);
% グラフの配置調整
figure('Visible', 'off', 'Position', [100, 100, 1800, 600], 'Color', 'white');
ax1 = axes('Position', [0.05, 0.1, 0.28, 0.8]);
ax2 = axes('Position', [0.37, 0.1, 0.28, 0.8]);
ax3 = axes('Position', [0.69, 0.1, 0.28, 0.8]);
[Y, X] = meshgrid(yl, xl); % X, Yグリッドを生成
% レンジの決定
xmin = min(xl);
xmax = max(xl);
ymin = min(yl);
ymax = max(yl);
% テキストボックスのテキスト作成
str_1 = ['$ \bar{a} = ', num2str(abar, '%.1e'), ...
', \, W = ', num2str(w, '%.1e'), ...
', \, M_p = ', num2str(mp, '%.1f'), '$'];
str_2 = ['$ \xi = ', num2str(xi, '%.2f'), ...
', \, k = ', num2str(k, '%.1f'), ...
', \, \theta_0 = ', num2str(theta_0, '%.1f'), ...
', \, \Delta = ', num2str(sc, '%.2f'), '$'];
str = {str_1, str_2};
% テキストボックスの位置とサイズを定義
dim1 = [0.12 0.1 0.18 0.12];
dim2 = [0.44 0.1 0.18 0.12];
dim3 = [0.76 0.1 0.18 0.12];
% テキストボックスを図に追加
annotation('textbox', dim1, 'String', str, 'FitBoxToText', 'on', ...
'BackgroundColor', 'none', 'EdgeColor', 'none', 'FontSize', 12, 'Interpreter', 'latex');
annotation('textbox', dim2, 'String', str, 'FitBoxToText', 'on', ...
'BackgroundColor', 'none', 'EdgeColor', 'none', 'FontSize', 12, 'Interpreter', 'latex');
annotation('textbox', dim3, 'String', str, 'FitBoxToText', 'on', ...
'BackgroundColor', 'none', 'EdgeColor', 'none', 'FontSize', 12, 'Interpreter', 'latex');
% 初期条件 ----------------------------------------------------------------
r0 = 2 * dx;
for j = 1 : ny
for i = 1 : nx
r = sqrt((i*dx)^2 +((ny-j)*dy)^2) - r0;% 原点に核生成
%r = sqrt( ((i-(nx/2))*dx)^2 +((j-(ny/2)) * dy)^2 ) - r0;% 中央に核生成
phi(i, j) = 0.5 * (1 - tanh(sqrt(2 * w) / (2 * abar) * r)); % フェーズフィールド変数の初期分布
T(i,j,1) = T_0 + phi(i,j) * (T_melt - T_0); % 温度の初期分布
end
end
% 時間発展 ----------------------------------------------------------------
for step = 1 : maxstep
step
% 界面異方性を考慮した勾配係数計算のループ
for i = 1 : nx
for j = 1 : ny
% 隣接節点
ip = i + 1;
im = i - 1;
jp = j + 1;
jm = j - 1;
% ノイマン条件
if ip > nx
ip = nx;
end
if im < 1
im = 1;
end
if jp > ny
jp = ny;
end
if jm < 1
jm = 1;
end
% 法線ベクトル計算
nvec_x(i,j) = -(phi(ip,j) - phi(im,j))/(2*dx);
nvec_y(i,j) = -(phi(i,jp) - phi(i,jm))/(2*dy);
% 法線ベクトルとX軸が成す角度の計算。基底ベクトルの正負に応じて割り当てる角度範囲を変える。
if nvec_x(i,j)>0 && nvec_y(i,j)>0
theta(i,j) = atan(nvec_y(i,j)/nvec_x(i,j));
elseif nvec_x(i,j)<0 && nvec_y(i,j)>0
theta(i,j) = atan(nvec_y(i,j)/nvec_x(i,j)) + pi;
elseif nvec_x(i,j)<0 && nvec_y(i,j)<0
theta(i,j) = atan(nvec_y(i,j)/nvec_x(i,j)) + pi;
elseif nvec_x(i,j)>0 && nvec_y(i,j)<0
theta(i,j) = atan(nvec_y(i,j)/nvec_x(i,j)) + 2*pi;
elseif nvec_x(i,j)>0 && nvec_y(i,j)==0
theta(i,j) = 0;
elseif nvec_x(i,j)==0 && nvec_y(i,j)>0
theta(i,j) = pi/2;
elseif nvec_x(i,j)<0 && nvec_y(i,j)==0
theta(i,j) = pi;
elseif nvec_x(i,j)==0 && nvec_y(i,j)<0
theta(i,j) = 3*pi/2;
end
% 勾配エネルギー係数a
a(i,j) = abar * (1 + xi * cos(k*(theta(i,j)-theta_0)));
% 1階微分
dphidx(i,j) = (phi(ip,j)-phi(im,j))/(2*dx);
dphidy(i,j) = (phi(i,jp)-phi(i,jm))/(2*dy);
% 9点ラプラシアン
lap_phi(i,j) = (2 * (phi(ip,j) + phi(im,j) + phi(i,jp) + phi(i,jm)) ...
+ phi(ip,jp) + phi(im,jm) + phi(im,jp) + phi(ip,jm) - 12 * phi(i,j)) / (3 * dx^2);
lap_temp(i,j) = (2 * (T(ip,j) + T(im,j) + T(i,jp) + T(i,jm))...
+ T(ip,jp) + T(im,jm) + T(im,jp) + T(ip,jm) - 12 * T(i,j)) / (3 * dx^2);
end
end
% フェーズフィールド変数、温度場の時間発展ループ
for i = 1 : nx
for j = 1 : ny
% 隣接節点
ip = i + 1;
im = i - 1;
jp = j + 1;
jm = j - 1;
% ノイマン条件
if ip > nx
ip = nx;
end
if im < 1
im = 1;
end
if jp > ny
jp = ny;
end
if jm < 1
jm = 1;
end
% aの微係数計算
A(i, j) = a(i,j)*abar*k*xi*sin(k*(theta(i,j)-theta_0))*dphidy(i,j);
B(i, j) = a(i,j)*abar*k*xi*sin(k*(theta(i,j)-theta_0))*dphidx(i,j);
dadx(i,j) = (a(ip,j)-a(im,j))/(2*dx);
dady(i,j) = (a(i,jp)-a(i,jm))/(2*dy);
% 化学的駆動力の増分の計算
dg(i,j) = latent * (T(i,j)-T_melt)/T_melt;
% フェーズフィールド変数の増分の計算
% エネルギー密度関数１
% dphi(i,j) =...
% -6*dg(i,j)*phi(i,j)*(1-phi(i,j)) - 2*w*phi(i,j)*(1-phi(i,j))*(1-2*phi(i,j))...
% +2*a(i,j,step)*(dadx(i,j)*dphidx(i,j)+dady(i,j)*dphidy(i,j))...
% +a(i,j,step)^(2)*lap_phi(i,j)...
% +(A(ip,j)-A(im,j))/(2*dx)-((B(i,jp)-B(i,jm))/(2*dy));
% エネルギー密度関数２
noize = 0.1*(2 * rand - 1);
dphi(i,j) =...
-30*dg(i,j)*phi(i,j)^2*(1-phi(i,j))^2 - 2*w*phi(i,j)*(1-phi(i,j))*(1-2*phi(i,j))...
+2*a(i,j)*(dadx(i,j)*dphidx(i,j)+dady(i,j)*dphidy(i,j))...
+a(i,j)^(2)*lap_phi(i,j)...
+(A(ip,j)-A(im,j))/(2*dx)-((B(i,jp)-B(i,jm))/(2*dy))...
+ 4 * w * phi(i,j)*(1-phi(i,j))* noize;
% フェーズフィールド変数の時間発展
phi_new(i,j) = phi(i,j) + dt * mp * dphi(i,j);
% 温度場の増分の計算
%q(i,j) = 6*latent*phi(i,j)*(1-phi(i,j))*(phi_new(i,j)-phi(i,j))/dt;
q(i,j) = 30*latent*phi(i,j)^2*(1-phi(i,j))^2*(phi_new(i,j)-phi(i,j))/dt;% 9点ラプラシアン
dT(i,j) = (1 / cv)*(K* lap_temp(i,j) +q(i,j));
% 温度場の時間発展
T_new(i,j) = T(i,j) + dt *dT(i,j);
end
end
% タイムステップを置き換える
phi = phi_new;
T = T_new;
% 動画保存
if step == 1
plotconfig(xl, yl, xmin,xmax, ymin,ymax,step,dt,ax1,ax2,ax3,T_0,T_melt,phi,T,q)
v = VideoWriter(filename,'MPEG-4');
v.FrameRate = 100;
open(v);
elseif mod(step, 60) == 0
plotconfig(xl, yl, xmin,xmax, ymin,ymax,step,dt,ax1,ax2,ax3,T_0,T_melt,phi,T,q)
frame = getframe(gcf);
writeVideo(v,frame);
end
end
% 動画ファイルを閉じる
close(v);
% 以下ローカル関数---------------------------------------------------------------
% 可視化
function [] = plotconfig(xl, yl, xmin,xmax, ymin,ymax,step,dt,ax1,ax2,ax3,T_0,T_melt,phi,T,q)
% Phase-field サブプロット
axes(ax1);
cla(ax1);
imagesc(xl, yl, imrotate(phi(:,:), 90), [0, 1]);
title(ax1, sprintf('phase-field $\\phi$ at %.3e [s]', step*dt), 'Interpreter', 'latex');
cb1 = colorbar(ax1);
title(cb1, '$\phi$', 'Interpreter', 'latex');
axis(ax1, 'equal');
axis(ax1, 'tight');
xlabel('$$x\rm{[m]}$$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('$$y\rm{[m]}$$', 'Interpreter', 'latex', 'FontSize', 14);
xticks(linspace(xmin, xmax, 6)); % x軸の目盛りを設定
yticks(linspace(ymin, ymax, 6)); % y軸の目盛りを設定
ax1.XAxis.TickLabelFormat = '%,.1f';
ax1.YAxis.TickLabelFormat = '%,.1f';
set(ax1, 'YDir', 'normal', 'FontName', 'Times New Roman', 'FontSize', 12,'LineWidth',1.2);
% 温度場 サブプロット
T_min = T_0;
T_max = T_melt;
axes(ax2);
cla(ax2);
imagesc(xl, yl, imrotate(T(:,:), 90), [T_min, T_max]);
title(ax2, sprintf('temperature $T \\rm{[K]}$ at %.3e [s]', step*dt), 'Interpreter', 'latex');
cb2 = colorbar(ax2);
title(cb2, '$T$', 'Interpreter', 'latex');
axis(ax2, 'equal');
axis(ax2, 'tight');
xlabel('$$x\rm{[m]}$$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('$$y\rm{[m]}$$', 'Interpreter', 'latex', 'FontSize', 14);
xticks(linspace(xmin, xmax, 6)); % x軸の目盛りを設定
yticks(linspace(ymin, ymax, 6)); % y軸の目盛りを設定
ax2.XAxis.TickLabelFormat = '%,.1f';
ax2.YAxis.TickLabelFormat = '%,.1f';
set(ax2, 'YDir', 'normal', 'FontName', 'Times New Roman', 'FontSize', 12,'LineWidth',1.2);
% 潜熱コンター
q_min = -3.106811707055531e+17;
q_max = 2.210685380238862e+19;
axes(ax3);
cla(ax3);
imagesc(xl, yl, imrotate(q(:,:), 90), [q_min, q_max]);
title(ax3, sprintf('latent heat $q \\rm{[W/m^3]}$ at %.3e [s]', step*dt), 'Interpreter', 'latex');
cb3 = colorbar(ax3);
title(cb3, '$q$', 'Interpreter', 'latex');
axis(ax3, 'equal');
axis(ax3, 'tight');
xlabel('$$x\rm{[m]}$$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('$$y\rm{[m]}$$', 'Interpreter', 'latex', 'FontSize', 14);
xticks(linspace(xmin, xmax, 6)); % x軸の目盛りを設定
yticks(linspace(ymin, ymax, 6)); % y軸の目盛りを設定
ax3.XAxis.TickLabelFormat = '%,.1f';
ax3.YAxis.TickLabelFormat = '%,.1f';
set(ax3, 'YDir', 'normal', 'FontName', 'Times New Roman', 'FontSize', 12, 'LineWidth', 1.2);
% 新しいプロットの時、軸設定を保持したまま前のグラフィックスオブジェクトを消去
set(gca,'nextplot','replacechildren');
end解析結果
以下に3条件の解析結果をまとめて示します。各動画は左から、フェーズフィールド変数、温度、潜熱のコンター図を示します。
過冷度が大きくなるにつれて、成長速度が大きくなり、アームの発達が活発化している様子が確認できます。これは過冷度の増加に伴い、化学的駆動力\(l\left(\frac{T-T_{m}}{T_{m}}\right)\)が増加していることを示しています。
また、\(\Delta = 0.3\)では、デンドライト構造が発生しませんでした。これは過冷度が小さいため、化学的駆動力\(l\left(\frac{T-T_{m}}{T_{m}}\right)\)が小さくなり、微小凹凸が発達しないためと考えられます。
おわりに
今回は純物質凝固モデルに化学的駆動力のゆらぎを考慮することで、デンドライト構造のシミュレーションを行いました。次回からは2元合金の凝固についてまとめていきたいと思います。
参考文献
[1] 高木 知弘ら, “フェーズフィールド法”, 養賢堂, 2012/3/2
[2] 山中 晃徳ら, “Pythonによるフェーズフィールド法入門 基礎理論からデータ同化の実装まで”, 丸善出版, 2023/12/15
[3] 小山 敏幸ら, “フェーズフィールド法入門 (計算力学レクチャーコース) “, 丸善出版, 2013/4/13
コメント