はじめに
前回記事では、1次元アレン-カーン方程式を解き、各パラメーターが解に与える影響について考察しました。一方で、デンドライド成長などで考慮される界面異方性は2次元以上で生じるエネルギー効果です。そこで今回は2次元アレンカーン方程式を実際に解き、前回同様、各パラメーターが解に与える影響についてまとめておきます。
2次元アレン-カーン方程式の離散化
汎関数微分を計算したアレン-カーン方程式は次のように書かれました。(過去の記事参照)
\frac{\partial \phi}{\partial t}&=M_{\phi}\left[a^2\nabla^2 \phi-\left(g_{a}-g_{b}\right)\frac{\partial p(\phi)}{\partial \phi}-W\frac{\partial q(\phi)}{\partial \phi}\right] \tag{1}\\
p(\phi)&=\phi ^2(3-2\phi) \tag{2}\\
q(\phi)&=\phi ^2 (1-\phi)^2 \tag{3} \\
\end{align}
前回記事と同様に式\((2),(3)\)を用いると、式\((1)\)右辺の第二、第三項は次のように計算できます。
\frac{\partial p(\phi)}{\partial \phi} &=6\phi\left(1-\phi\right) \tag{4}\\
\end{align}
\frac{\partial q(\phi)}{\partial \phi} &= 2\phi (1-\phi)(1-2\phi) \tag{5}\\
\end{align}
これらを式\((1)\)に代入し、二次元として考えると、以下の\(\phi\)に関する方程式を得ることができます。
\frac{\partial \phi}{\partial t}&=M_{\phi}\left[a^2\left(\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}\right) -6\left(g_{a}-g_{b}\right)\phi\left(1-\phi\right)-2W\phi (1-\phi)(1-2\phi) \right] \tag{6}\\
\end{align}
さらに式\((6)\)を離散化していきます。左辺の時間微分はオイラー陽解法で、右辺の空間2階微分は中心差分法で離散化すると次式が得られます。
&\frac{\phi^{n+1}_{i,j}-\phi^n_{i,j}}{\Delta t}=M_{\phi}\left[a^2\left(\frac{\phi^{n}_{i-1,j}-2\phi^{n}_{i,j}+\phi^{n}_{i+1,j}}{\left(\Delta x\right)^2}+\frac{\phi^{n}_{i,j-1}-2\phi^{n}_{i,j}+\phi^{n}_{i,j+1}}{\left(\Delta y\right)^2}\right)\\
\qquad \qquad \qquad \qquad -6\left(g_{a}-g_{b}\right)\phi ^{n} _{i,j}\left(1-\phi^{n} _{i,j}\right)-2W\phi_{i,j}^{n} (1-\phi_{i,j}^{n})(1-2\phi_{i,j}^{n}) \right] \\
&\rightarrow
\phi^{n+1}_{i}=\phi^n_{i}\\
&\quad \quad \quad+{\Delta t}M_{\phi}\left[a^2\left(\frac{\phi^{n}_{i-1,j}-2\phi^{n}_{i,j}+\phi^{n}_{i+1,j}}{\left(\Delta x\right)^2}+\frac{\phi^{n}_{i,j-1}-2\phi^{n}_{i,j}+\phi^{n}_{i,j+1}}{\left(\Delta y\right)^2}\right)\\
\qquad \qquad \qquad \qquad -6\Delta g\phi ^{n} _{i}\left(1-\phi^{n} _{i}\right)-2W\phi^{n}_{i} (1-\phi^{n}_{i})(1-2\phi^{n}_{i}) \right] \tag{7}\\
\end{align}
ここで\(\Delta g =g_{a}-g_{b}\)としました。以上が2次元アレン-カーン方程式の離散化形式です。
初期条件・境界条件
初期条件として、原点を中心に半径1の円領域を境に相a(\(\phi=1\))から相b(\(\phi=0\))に連続的に変化する以下の関数を考えます。
\phi _0&=\frac{1}{2}\left(1-\tanh(2r)\right)\tag{8}\\
r &= \left(x^2 + y^2\right)^{\frac{1}{2}} – 1\tag{9}\\
\end{align}

また、境界条件は両端で勾配が0になるノイマン条件を設定することとします。
MATLABコード
以下に今回使用したMATLABコードを記載します。
% 初期化処理 ---------------------------------------------------
clc; close all;
clear all;
% パラメーター設定 ---------------------------------------------
dx = 0.1; % 差分格子点の間隔[m]
dy = 0.1; % 差分格子点の間隔[m]
a = 1; % 勾配エネルギ―係数[(J/m)^1/2]
w = 0; % エネルギー障壁[J/m3]
dg = 0; % 化学的駆動力[J/m3]
mp = 1; % 界面モビリティ[m3/js]
dt = 0.001; % タイムステップ[s]
maxstep = 5000; % 最大ステップ数
xmax = 3;
xmin = -3;
ymax = 3;
ymin = -3;
xl = xmin:dx:xmax;
yl = ymin:dy:ymax;
% 変数宣言
nx = size(xl,2);% 差分格子点数
ny = size(yl,2);% 差分格子点数
phi = zeros(nx,ny,maxstep);
dphi = zeros(nx,ny,maxstep);
% 安定条件の確認(cが0.5以下でないと不安定)
c = 2*a^2 *dt / (dx^2+dy^2);
% 初期条件 ---------------------------------------------------
r_nuclei = 1;% 原点から半径r0=r_nucleiの領域をφ=1とする。
for i = 1 : nx
for j = 1 : ny
r = sqrt(xl(i)^2 + yl(j)^2) - r_nuclei;
phi(i,j,1) = 0.5 * (1 - tanh(2*r)); % フェーズフィールド変数の初期分布の設定
end
end
% 時間発展 ---------------------------------------------------
for step = 1 : maxstep - 1
for i = 2 : nx - 1
for j = 2 : ny - 1
dphi(i,j,step) = a^2 * ((phi(i-1,j,step) -2*phi(i,j,step) +phi(i+1,j,step))/dx^2 ...
+(phi(i,j-1,step) -2*phi(i,j,step) +phi(i,j+1,step))/dy^2) ...
- dg * 6 *phi(i,j,step)*(1-phi(i,j,step)) - 2*w*phi(i,j,step)*(1-phi(i,j,step))*(1-2*phi(i,j,step));
phi(i,j,step + 1) = phi(i,j,step) + dt * mp * dphi(i,j,step);
end
end
% ノイマン境界条件
phi(1,:,step+1) = phi(2,:,step+1);
phi(nx,:,step+1) = phi(nx-1,:,step+1);
phi(:,1,step+1) = phi(:,2,step+1);
phi(:,ny,step+1) = phi(:,ny-1,step+1);
end
% 可視化 -----------------------------------------------------
% VideoWriter オブジェクトの初期化
filename = ['a = ', num2str(a, '%.0f'), ', w = ', num2str(w, '%.0f'), ...
', dg = ', num2str(dg, '%.0f'), ', mp = ', num2str(mp, '%.3f'),'.mp4'];
outputVideo = VideoWriter(filename, 'MPEG-4');
outputVideo.FrameRate = 30; % フレームレートの設定
open(outputVideo);
% グラフィックオブジェクトの初期化と再利用
fig = figure(1);
plotHandle = imagesc(xl, yl, phi(:,:,1)); % 初期プロット
% 軸ラベルとタイトルの設定
xlabel('$$x$$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('$$y$$', 'Interpreter', 'latex', 'FontSize', 14);
% Y軸の向きを正に設定
set(gca, 'YDir', 'normal');
% 軸の設定
ax = gca;
ax.FontName = 'Times New Roman';
ax.FontSize = 12;
ax.LineWidth = 1.2; % 枠の太さ
axis equal;
axis tight;
axis on;
% カラーバーを追加
cb = colorbar;
title(cb, ' $\phi$', 'Interpreter', 'latex');
caxis([0 1])
% 背景を白に
fig=gcf;
fig.Color='white';
% タイトルと凡例の設定
titleHandle = title(sprintf('time = %.3f [s]', 0));
% テキストボックスの設定
str = ['$a = ', num2str(a, '%.1f'), '\, , W = ', num2str(w, '%.1f'), ', \, \Delta g = ', num2str(dg, '%.1f'), '$'];
dim = [0.43 0.1 0.2 0.1];% テキストボックスの位置とサイズ [左下のX座標, 左下のY座標, 幅, 高さ]
annotation('textbox', dim, 'String', str, 'FitBoxToText', 'on', 'BackgroundColor', ...
'white', 'FontSize', 10, 'Interpreter', 'latex');% テキストボックスを図に追加
% 動画のフレームを記録するためのループ
for step = 1 : maxstep
if mod(step, 10) == 0
% データ点の更新。imagescで作成されたイメージオブジェクトの'CData'を更新します。
set(plotHandle, 'CData', phi(:,:,step));
% タイトルの更新
set(titleHandle, 'String', sprintf('time = %.3f [s]', step * dt));
drawnow; % グラフを更新
frame = getframe(gcf); % 現在のフレームをキャプチャ
writeVideo(outputVideo, frame); % ビデオファイルに書き込み
end
end
% 動画の書き出しを完了
close(outputVideo);
各パラメーターの影響
前回記事同様、勾配係数\(a\)、エネルギー障壁\(W\)、化学的駆動力\(\Delta g\)が解にどのような影響を与えるか確認します。
まず、勾配係数\(a\)の影響を見るために、以下のようなパラメーターで計算を実行します。
パラメーター | 値 |
---|---|
\(a\) |      1      |
\(W\) | 0 |
\(\Delta g\) | 0 |
1次元の時と同様、時間とともに解が鈍っていき、やがて消失します。これは拡散方程式の解と同等です。
次にエネルギー障壁\(W\)の影響を見るために、以下のようなパラメーターで計算を実行します。
パラメーター | 値 |
---|---|
\(a\) | 0 |
\(W\) |      1      |
\(\Delta g\) | 0 |
初期状態で緩やかだった解が時間とともにステップ型に変化します。これはエネルギーがより低い相状態に解を移行させる、ダブルウェルポテンシャルの特徴です。
最後に、化学的駆動力\(\Delta g=g_{a}-g_{b}\)の影響を見るために、次のパラメーターで計算を実行します。
パラメーター | 値 |
---|---|
\(a\) | 0 |
\(W\) | 0 |
\(\Delta g\) |      1      |
\(\Delta g=g_{a}-g_{b}=1>0\)より、時間とともによりエネルギーが大きいa相(\(\phi =1\))が縮小し、b相(\(\phi =0\))に置き換わっていきます。これも1次元の時と同様の特徴です。
おわりに
今回は2次元アレン-カーン方程式を解き、各項が与える影響についてまとめました。次回からは界面異方性について考えていきたいと思います。
参考文献
[1] 高木 知弘ら, “フェーズフィールド法”, 養賢堂, 2012/3/2
[2] 山中 晃徳ら, “Pythonによるフェーズフィールド法入門 基礎理論からデータ同化の実装まで”, 丸善出版, 2023/12/15
コメント