はじめに
前回は以下の界面異方性を考慮した2次元アレン-カーン方程式を導出しました。
\begin{align}
\frac{\partial \phi}{\partial t}&=-M_{\phi}\left[6\left(g_{a}-g_{b}\right)\phi\left(1-\phi\right)+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{1}\\
a(\theta)&=\bar{a}\left[1+\xi \cos \left\{ k\left( \theta -\theta _0\right)\right\} \right]
\tag{2}
\end{align}
\frac{\partial \phi}{\partial t}&=-M_{\phi}\left[6\left(g_{a}-g_{b}\right)\phi\left(1-\phi\right)+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{1}\\
a(\theta)&=\bar{a}\left[1+\xi \cos \left\{ k\left( \theta -\theta _0\right)\right\} \right]
\tag{2}
\end{align}
今回は上式を離散化し、シミュレーションしてみたいと思います。
離散化
式\((1)\)を離散化していきます。左辺の時間微分はオイラー陽解法、右辺の空間微分は中心差分法を利用します。
\begin{align}
\frac{\phi_{i,j}^{n+1}-\phi_{i,j}^{n}}{\Delta t}&=-M_{\phi}\left[-6\left(g_{a}-g_{b}\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. \quad\quad\quad\quad\quad-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\left(\frac{\phi_{i+1,j}^{n}-2\phi_{i,j}^{n}+\phi_{i-1,j}^{n}}{\left(\Delta x\right)^2}+\frac{\phi_{i,j+1}^{n}-2\phi_{i,j}^{n}+\phi_{i,j-1}^{n}}{\left(\Delta y\right)^2} \right)
\right.\\
&\left.\quad\quad\quad\quad\quad\quad\quad\quad
-\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}
\right] \\
\rightarrow
\phi_{i,j}^{n+1}&=\phi_{i,j}^{n}+M_{\phi}\Delta t\left[-6\left(g_{a}-g_{b}\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. \quad\quad\quad\quad\quad+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\left(\frac{\phi_{i+1,j}^{n}-2\phi_{i,j}^{n}+\phi_{i-1,j}^{n}}{\left(\Delta x\right)^2}+\frac{\phi_{i,j+1}^{n}-2\phi_{i,j}^{n}+\phi_{i,j-1}^{n}}{\left(\Delta y\right)^2} \right)
\right.\\
&\left.\quad\quad\quad\quad\quad\quad\quad\quad
+\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}
\right] \tag{3}\\
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{4}\\
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{5}\\
a_{i,j}^{n}&=\bar{a}\left[1+\xi \cos \left\{ k\left( \theta_{i,j}^{n} -\theta _0\right)\right\} \right]
\tag{6}
\end{align}
\frac{\phi_{i,j}^{n+1}-\phi_{i,j}^{n}}{\Delta t}&=-M_{\phi}\left[-6\left(g_{a}-g_{b}\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. \quad\quad\quad\quad\quad-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\left(\frac{\phi_{i+1,j}^{n}-2\phi_{i,j}^{n}+\phi_{i-1,j}^{n}}{\left(\Delta x\right)^2}+\frac{\phi_{i,j+1}^{n}-2\phi_{i,j}^{n}+\phi_{i,j-1}^{n}}{\left(\Delta y\right)^2} \right)
\right.\\
&\left.\quad\quad\quad\quad\quad\quad\quad\quad
-\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}
\right] \\
\rightarrow
\phi_{i,j}^{n+1}&=\phi_{i,j}^{n}+M_{\phi}\Delta t\left[-6\left(g_{a}-g_{b}\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. \quad\quad\quad\quad\quad+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\left(\frac{\phi_{i+1,j}^{n}-2\phi_{i,j}^{n}+\phi_{i-1,j}^{n}}{\left(\Delta x\right)^2}+\frac{\phi_{i,j+1}^{n}-2\phi_{i,j}^{n}+\phi_{i,j-1}^{n}}{\left(\Delta y\right)^2} \right)
\right.\\
&\left.\quad\quad\quad\quad\quad\quad\quad\quad
+\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}
\right] \tag{3}\\
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{4}\\
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{5}\\
a_{i,j}^{n}&=\bar{a}\left[1+\xi \cos \left\{ k\left( \theta_{i,j}^{n} -\theta _0\right)\right\} \right]
\tag{6}
\end{align}
ここで\(\theta_{i,j}^{n}\)は時間・空間分布を有し、次式で与えられます。法線ベクトル\(\boldsymbol{n}_{i,j}^{n}\)の向きに応じて計算式が変化することに注意しましょう。
\begin{align}
\theta_{i,j}^{n}&=
\begin{cases}
\tan ^{-1}\left(\frac{ny_{i,j}^{n}}{nx_{i,j}^{n}}\right) & nx_{i,j}^{n}> 0かつny_{i,j}^{n}> 0の場合\\
\tan ^{-1}\left(\frac{ny_{i,j}^{n}}{nx_{i,j}^{n}}\right) +\pi & nx_{i,j}^{n}< 0かつny_{i,j}^{n}> 0の場合\\
\tan ^{-1}\left(\frac{ny_{i,j}^{n}}{nx_{i,j}^{n}}\right) +\pi & nx_{i,j}^{n}< 0かつny_{i,j}^{n}< 0の場合\\
\tan ^{-1}\left(\frac{ny_{i,j}^{n}}{nx_{i,j}^{n}}\right) +2\pi & nx_{i,j}^{n}> 0かつny_{i,j}^{n}< 0の場合\\
0 & nx_{i,j}^{n}> 0かつny_{i,j}^{n}= 0の場合\\
\frac{\pi}{2} & nx_{i,j}^{n}= 0かつny_{i,j}^{n}> 0の場合\\
\pi & nx_{i,j}^{n}< 0かつny_{i,j}^{n}= 0の場合\\
\frac{3\pi}{2} & nx_{i,j}^{n}=0かつny_{i,j}^{n}< 0の場合\\
\end{cases}
\tag{7}\\
\boldsymbol{n}_{i,j}^{n}
&=
\begin{pmatrix}
nx_{i,j}^{n} \\ ny_{i,j}^{n}
\end{pmatrix}
= -\nabla \phi
= –
\begin{pmatrix}
\frac{\phi_{i+1,j}^{n}-\phi_{i-1,j}^{n}}{2\Delta x} \\ \frac{\phi_{i,j+1}^{n}-\phi_{i,j-1}^{n}}{2\Delta y}
\end{pmatrix}
\tag{8}
\end{align}
\theta_{i,j}^{n}&=
\begin{cases}
\tan ^{-1}\left(\frac{ny_{i,j}^{n}}{nx_{i,j}^{n}}\right) & nx_{i,j}^{n}> 0かつny_{i,j}^{n}> 0の場合\\
\tan ^{-1}\left(\frac{ny_{i,j}^{n}}{nx_{i,j}^{n}}\right) +\pi & nx_{i,j}^{n}< 0かつny_{i,j}^{n}> 0の場合\\
\tan ^{-1}\left(\frac{ny_{i,j}^{n}}{nx_{i,j}^{n}}\right) +\pi & nx_{i,j}^{n}< 0かつny_{i,j}^{n}< 0の場合\\
\tan ^{-1}\left(\frac{ny_{i,j}^{n}}{nx_{i,j}^{n}}\right) +2\pi & nx_{i,j}^{n}> 0かつny_{i,j}^{n}< 0の場合\\
0 & nx_{i,j}^{n}> 0かつny_{i,j}^{n}= 0の場合\\
\frac{\pi}{2} & nx_{i,j}^{n}= 0かつny_{i,j}^{n}> 0の場合\\
\pi & nx_{i,j}^{n}< 0かつny_{i,j}^{n}= 0の場合\\
\frac{3\pi}{2} & nx_{i,j}^{n}=0かつny_{i,j}^{n}< 0の場合\\
\end{cases}
\tag{7}\\
\boldsymbol{n}_{i,j}^{n}
&=
\begin{pmatrix}
nx_{i,j}^{n} \\ ny_{i,j}^{n}
\end{pmatrix}
= -\nabla \phi
= –
\begin{pmatrix}
\frac{\phi_{i+1,j}^{n}-\phi_{i-1,j}^{n}}{2\Delta x} \\ \frac{\phi_{i,j+1}^{n}-\phi_{i,j-1}^{n}}{2\Delta y}
\end{pmatrix}
\tag{8}
\end{align}
以上が界面異方性を考慮したアレン-カーン方程式の離散化形式です。
計算フローは次の通りです。
以上のように界面異方性を考慮すると、\(\phi\)の計算結果に応じて勾配係数\(a(\theta)\)(=拡散しやすい方向)が常に変化するため、異方性を考慮しない時と比べて解が飛躍的に複雑になります。
解の安定性条件
アレン-カーン方程式の解の安定性条件はどのように考えるべきでしょうか？ここでは拡散方程式との比較から簡易的に考慮することとします。
一般に、拡散係数\(D\)で表される2次元拡散方程式は次式で与えられました。
\begin{align}
\frac{\partial \phi}{\partial t}&=D\left(\frac{\partial ^2\phi}{\partial x^2}+\frac{\partial ^2\phi}{\partial y^2} \right)\tag{9}
\end{align}
\frac{\partial \phi}{\partial t}&=D\left(\frac{\partial ^2\phi}{\partial x^2}+\frac{\partial ^2\phi}{\partial y^2} \right)\tag{9}
\end{align}
上式を陽解法で解く場合の解の安定性条件は次のように表されます。
\begin{align}
\Delta t\leq \frac{1}{2D}\frac{\left(\Delta x\right)^2\left(\Delta y\right)^2}{\left(\Delta x\right)^2+\left(\Delta y\right)^2}\tag{10}
\end{align}
\Delta t\leq \frac{1}{2D}\frac{\left(\Delta x\right)^2\left(\Delta y\right)^2}{\left(\Delta x\right)^2+\left(\Delta y\right)^2}\tag{10}
\end{align}
一方で、アレン-カーン方程式\((1)\)は非線形項を含みますが、\(M_{\phi}a^2\)を拡散係数とする拡散方程式とみなすこともできます。したがって、拡散項の影響だけ考慮した場合、解の安定性条件は次式で与えてもよいでしょう。
\begin{align}
\Delta t\leq \frac{1}{2M_{\phi}a^2}\frac{\left(\Delta x\right)^2\left(\Delta y\right)^2}{\left(\Delta x\right)^2+\left(\Delta y\right)^2}\tag{11}
\end{align}
\Delta t\leq \frac{1}{2M_{\phi}a^2}\frac{\left(\Delta x\right)^2\left(\Delta y\right)^2}{\left(\Delta x\right)^2+\left(\Delta y\right)^2}\tag{11}
\end{align}
ここで勾配係数\(a\)は界面異方性を持つため、法線ベクトルの向き\(\theta\)によって変化します。この為、安定性条件としては\(a\)の最大値\(\max(a)\)を用いればよいでしょう。
\begin{align}
\Delta t\leq \frac{1}{2M_{\phi}\cdot\max(a)^2}\frac{\left(\Delta x\right)^2\left(\Delta y\right)^2}{\left(\Delta x\right)^2+\left(\Delta y\right)^2}\tag{12}
\end{align}
\Delta t\leq \frac{1}{2M_{\phi}\cdot\max(a)^2}\frac{\left(\Delta x\right)^2\left(\Delta y\right)^2}{\left(\Delta x\right)^2+\left(\Delta y\right)^2}\tag{12}
\end{align}
無次元数の形にすると
\begin{align}
c= 2M_{\phi}\Delta t\cdot\max(a)^2 \frac{\left(\Delta x\right)^2+\left(\Delta y\right)^2}{\left(\Delta x\right)^2\left(\Delta y\right)^2}\leq 1\tag{13}
\end{align}
c= 2M_{\phi}\Delta t\cdot\max(a)^2 \frac{\left(\Delta x\right)^2+\left(\Delta y\right)^2}{\left(\Delta x\right)^2\left(\Delta y\right)^2}\leq 1\tag{13}
\end{align}
と書かれ、\(c\)が1以下になるような\(\Delta t\)が安定条件です。式\((2)\)より\(a\)の最大値を大きくする\(\bar{a},\xi\)といったパラメーターの値には特に注意する必要があります。
以上がアレン-カーン方程式の簡略化された解の安定性条件です。より厳密な条件を求めるには拡散項以外の非線形項も考慮し、安定性解析を行う必要があると考えられます。しかしながら、そのような計算は本記事の本筋ではないので今回は上記のように簡単化して扱うこととします。実際には式\((13\))に\(\frac{1}{2}\)を掛けておくなど、安定側の係数をかけて用いるのが良いでしょう。
初期条件・境界条件
初期条件として、原点を中心に半径１の円領域を境に相a（\(\phi=1\)）から相b（\(\phi=0\)）に連続的に変化する以下の関数を考えます。
\begin{align}
\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}
\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}
また境界条件は両端で勾配が０になるノイマン条件を設定することとします。
計算結果
では実際に計算してみましょう。今回は異方性強度\(\xi\)、異方性モード\(k\)を変化させた場合、フェーズフィールド変数\(\phi\)、勾配エネルギー係数\(a\)、界面法線ベクトル\(n\)がどのように変化するか確認します。
異方性とは関係ないパラメーター\(M_{\phi},W,\Delta g=g_{a}-g_{b}\)は\(M_{\phi}=1,W=0,\Delta g=0\)で固定しておくこととします。また、式\((2)\)より\(\bar{a}\)は勾配エネルギー係数をスケール倍する効果であり、比較的影響が分かりやすいので\(\bar{a}=0.5\)で固定しておきます。
まず、異方性強度\(\xi\)について\(\xi＝0.5,1.0,1.5\)と変えた場合、以下の結果が得られます。
以上より、\(\xi\)を大きくすると、勾配エネルギー係数\(a\)が\(x\)方向に集中し、\(x\)方向で\(\phi\)の拡散が強くなる様子が確認できます。したがって、異方性強度\(\xi\)は特定方向の異方性を強める働きがあることが分かります。
次に、異方性モード\(k\)について\(k＝1,2,3\)と変えた場合、以下の結果が得られます。
以上より、\(k\)を大きくすると、勾配エネルギー係数\(a\)の大きい領域が\(\frac{360}{k}\rm{[deg]}\)方向ごとに現れ、その方向で\(\phi\)の拡散が強くなる様子が確認できます。したがって、異方性モード\(k\)は異方性の周波数を高くし、次数を上げる働きがあることが分かります。
MATLABコード
以下に今回使用したMATLABコードを記載します。
% 初期化処理 --------------------------------------------------------------
clc; close all;
clear all;
% パラメーター設定 --------------------------------------------------------
% 物性設定
w = 0; % エネルギー障壁[J/m3]
dg = 0; % 化学的駆動力[J/m3]
mp = 1; % 界面モビリティ[m3/js]
abar = 0.5; % 基準となる勾配係数
xi = 0.5; % 異方性強度
k = 1; % 異方性モード
theta_0 = 0; % 優先成長方向
% 時間・空間ステップ設定
dt = 0.001; % タイムステップ[s]
maxstep = 3000; % 最大ステップ数
dx = 0.1; % 差分格子点の間隔[m]
dy = 0.1; % 差分格子点の間隔[m]
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);
nvec_x = zeros(nx,ny,maxstep);
nvec_y = zeros(nx,ny,maxstep);
theta = zeros(nx,ny);
thetadeg = zeros(nx,ny);
a = zeros(nx,ny,maxstep);
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);
% 安定条件の確認
theta_test = 0:2*pi/100:2*pi;
atest = max(abar * (1 + xi * cos(k*(theta_test-theta_0))));% 勾配係数の最大値を計算する。
c = (1/(2*atest^2*mp))*(dx^2*dy^2)/(dx^2+dy^2);
if dt > c
error('タイムステップが安定条件を満たしていません。')
end
% 初期条件 ----------------------------------------------------------------
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
step
% 界面異方性を考慮した勾配係数計算のループ
for i = 2 : nx - 1
for j = 2 : ny - 1
% 法線ベクトル計算
nvec_x(i,j,step) = -(phi(i+1,j,step) - phi(i-1,j,step))/(2*dx);
nvec_y(i,j,step) = -(phi(i,j+1,step) - phi(i,j-1,step))/(2*dy);
% 法線ベクトルとX軸が成す角度の計算。基底ベクトルの正負に応じて割り当てる角度範囲を変える。
if nvec_x(i,j,step)>0 && nvec_y(i,j,step)>0
theta(i,j) = atan(nvec_y(i,j,step)/nvec_x(i,j,step));
elseif nvec_x(i,j,step)<0 && nvec_y(i,j,step)>0
theta(i,j) = atan(nvec_y(i,j,step)/nvec_x(i,j,step)) + pi;
elseif nvec_x(i,j,step)<0 && nvec_y(i,j,step)<0
theta(i,j) = atan(nvec_y(i,j,step)/nvec_x(i,j,step)) + pi;
elseif nvec_x(i,j,step)>0 && nvec_y(i,j,step)<0
theta(i,j) = atan(nvec_y(i,j,step)/nvec_x(i,j,step)) + 2*pi;
elseif nvec_x(i,j,step)>0 && nvec_y(i,j,step)==0
theta(i,j) = 0;
elseif nvec_x(i,j,step)==0 && nvec_y(i,j,step)>0
theta(i,j) = pi/2;
elseif nvec_x(i,j,step)<0 && nvec_y(i,j,step)==0
theta(i,j) = pi;
elseif nvec_x(i,j,step)==0 && nvec_y(i,j,step)<0
theta(i,j) = 3*pi/2;
end
% 勾配エネルギー係数a
a(i,j,step) = abar * (1 + xi * cos(k*(theta(i,j)-theta_0)));
% 各種微係数の計算
dphidx(i,j) = (phi(i+1,j,step)-phi(i-1,j,step))/(2*dx);
dphidy(i,j) = (phi(i,j+1,step)-phi(i,j-1,step))/(2*dy);
lap_phi(i,j) = (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;
end
end
% フェーズフィールド変数の時間発展ループ
for i = 2 : nx - 1
for j = 2 : ny - 1
% aの微係数計算
A(i, j) = a(i,j,step)*abar*k*xi*sin(k*(theta(i,j)-theta_0))*dphidy(i,j);
B(i, j) = a(i,j,step)*abar*k*xi*sin(k*(theta(i,j)-theta_0))*dphidx(i,j);
dadx(i,j) = (a(i+1,j,step)-a(i-1,j,step))/(2*dx);
dady(i,j) = (a(i,j+1,step)-a(i,j-1,step))/(2*dy);
% フェーズフィールド変数の増分
dphi(i,j,step) =...
-6*dg*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))...
+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(i+1,j)-A(i-1,j))/(2*dx)-((B(i,j+1)-B(i,j-1))/(2*dy));
% 時間発展
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
% 可視化 ------------------------------------------------------------------
% ファイル名
paramstr =['abar = ', num2str(abar, '%.2f'), ', w = ', num2str(w, '%.2f'), ...
', dg = ', num2str(dg, '%.2f'), ', mp = ', num2str(mp, '%.2f'),...
', xi = ', num2str(xi, '%.2f'), ', k = ', num2str(k, '%.2f'),...
', theta_0 = ', num2str(theta_0, '%.2f')];
filename = [paramstr,'.mp4'];
% VideoWriter オブジェクトの初期化
outputVideo = VideoWriter(filename, 'MPEG-4');
outputVideo.FrameRate = 30; % フレームレートの設定
open(outputVideo);
% グラフの配置調整
figure(1);
set(gcf, 'Position', [100, 100, 1800, 600], 'Color', 'white');
ax1 = axes('Position', [0.03, 0.1, 0.28, 0.8]);
ax2 = axes('Position', [0.34, 0.1, 0.28, 0.8]);
ax3 = axes('Position', [0.65, 0.1, 0.28, 0.8]);
[Y, X] = meshgrid(yl, xl); % X, Yグリッドを生成
a_max = max(max(max(a)));
a_min = min(min(min(a)));
% テキストボックスのテキスト作成
str_1 = ['$ \bar{a} = ', num2str(abar, '%.1f'), ...
', \, W = ', num2str(w, '%.1f'), ...
', \, \Delta g = ', num2str(dg, '%.1f'), ...
', \, M_p = ', num2str(mp, '%.1f'), '$'];
str_2 = ['$ \xi = ', num2str(xi, '%.1f'), ...
', \, k = ', num2str(k, '%.1f'), ...
', \, \theta_0 = ', num2str(theta_0, '%.1f'), '$'];
str = {str_1, str_2};
% テキストボックスの位置とサイズを定義
dim1 = [0.11 0.1 0.18 0.12];
dim2 = [0.42 0.1 0.18 0.12];
dim3 = [0.76 0.07 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');
for nstep = 1 : maxstep - 1
nstep
if mod(nstep, 10) == 0
% Phase-field サブプロット
axes(ax1);
cla(ax1);
imagesc(xl, yl, imrotate(phi(:,:,nstep), 90), [0, 1]);
title(ax1, sprintf('phase-field $\\phi$ at %.3f [s]', nstep*dt), 'Interpreter', 'latex');
cb1 = colorbar(ax1);
title(cb1, '$\phi$', 'Interpreter', 'latex');
axis(ax1, 'equal');
axis(ax1, 'tight');
xlabel('$$x$$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('$$y$$', 'Interpreter', 'latex', 'FontSize', 14);
set(ax1, 'YDir', 'normal', 'FontName', 'Times New Roman', 'FontSize', 12,'LineWidth',1.2);
% 勾配エネルギー係数 サブプロット
axes(ax2);
cla(ax2);
imagesc(xl, yl, imrotate(a(:,:,nstep), 90), [a_min, a_max]);
title(ax2, sprintf('gradient energy coefficient $a$ at %.3f [s]', nstep*dt), 'Interpreter', 'latex');
cb2 = colorbar(ax2);
title(cb2, '$a$', 'Interpreter', 'latex');
axis(ax2, 'equal');
axis(ax2, 'tight');
xlabel('$$x$$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('$$y$$', 'Interpreter', 'latex', 'FontSize', 14);
set(ax2, 'YDir', 'normal', 'FontName', 'Times New Roman', 'FontSize', 12,'LineWidth',1.2);
% 法線ベクトル サブプロット
axes(ax3);
cla(ax3);
quiver(ax3, X, Y, nvec_x(:,:,nstep), nvec_y(:,:,nstep), 1.5);
title(ax3, sprintf('normal vector at %.3f [s]', nstep*dt));
axis(ax3, 'equal');
axis(ax3, 'tight');
xlabel('$$x$$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('$$y$$', 'Interpreter', 'latex', 'FontSize', 14);
set(ax3, 'YDir', 'normal', 'FontName', 'Times New Roman', 'FontSize', 12, 'LineWidth', 1.2);
% 画面を更新する
frame = getframe(gcf);
writeVideo(outputVideo, frame);
end
end
% 動画の書き出しを完了
close(outputVideo);
おわりに
今回は界面異方性を考慮したアレン-カーン方程式を実際に解き、パラメーターが与える影響について確認しました。
本記事では、簡単のためエネルギー障壁\(W\)、化学的駆動力\(\Delta g\)などのパラメーターは0として固定し、その影響を無視しました。当然ですがこれらの影響も考慮すると解の挙動はさらに複雑になります。しかしながら、これらのパラメーターが与える影響はこれまでに確認した通りなので、なんとなくどんな解になるかはイメージがつきやすいと思います。もし気になる場合は、実際にパラメーターを変えて計算を行ってみて下さい。
次回は純物質の凝固についてまとめたいと思います。
参考文献
[1] 高木 知弘ら, “フェーズフィールド法”, 養賢堂, 2012/3/2
[2] 山中 晃徳ら, “Pythonによるフェーズフィールド法入門 基礎理論からデータ同化の実装まで”, 丸善出版, 2023/12/15
コメント