Multi Phase Field Model with Anisotropic Interface Energy

Multi Phase Field Model with Anisotropic Interface Energy

Seunghwa Ryu, Yanming Wang and Wei Cai

This tutorial describes the theory and numerical implementation of a 3D multi phase phase field model with anisotropic interface energy. The model has been implemented in both matlab and in MD++. We describe the algorithm of these codes and how to use them in detail. A brief introduction of the phase field model can be found in Chapter 11 of

Computer Simulations of Dislocations, V. V. Bulatov and W. Cai, Oxford University Press (2006).

---

Download files

The matlab program can be downloaded from this link: Pf3d multi aniso.m. For the test case you need only one matlab file. You can download it into a directory of your choice, using the following command

wget http://micro.stanford.edu/mediawiki/images/f/fd/Pf3d_multi_aniso.m.txt -O Pf3d_multi_aniso.m

and run it within Matlab.

Theory

Similar to the case of multi phase field model with isotropic interface energy, the fundamental degrees of freedom and the corresponding free energy functional are in the same style. If we consider a system containing vapor, liquid and solid three phases, the only difference is that the gradient term in the free energy density functions of liquid-solid pair and solid-vapor pair now depends on the local orientation of the gradient of solid phase.

$F[\phi _{1}(\mathbf {x} ),...,\phi _{N}(\mathbf {x} )]=\int \sum _{k\neq {i}}{f_{ik}(\phi _{i}(\mathbf {x} ),\nabla \phi _{i}(\mathbf {x} ),\phi _{k}(\mathbf {x} ),\nabla \phi _{k}(\mathbf {x} ))d^{3}\mathbf {x} }$
Replacing $\varepsilon _{LS}$ with ${\frac {1}{2}}(\epsilon _{LS}(\mathbf {n_{S}} ))^{2}$, $\varepsilon _{SV}$ with ${\frac {1}{2}}(\epsilon _{SV}(\mathbf {n_{S}} ))^{2}$ and keeping $\varepsilon _{VL}$ a constant,
$f_{SV}(\phi _{S},\nabla \phi _{S},\phi _{V},\nabla \phi _{V})={\frac {1}{2}}(\epsilon _{SV}(\mathbf {n_{S}} ))^{2}(\phi _{V}\nabla \phi _{S}-\phi _{S}\nabla \phi _{V})^{2}+U_{SV}\phi _{S}^{2}\phi _{V}^{2}$
$f_{LS}(\phi _{L},\nabla \phi _{L},\phi _{S},\nabla \phi _{S})={\frac {1}{2}}(\epsilon _{LS}(\mathbf {n_{S}} ))^{2}(\phi _{S}\nabla \phi _{L}-\phi _{L}\nabla \phi _{S})^{2}+U_{LS}\phi _{L}^{2}\phi _{S}^{2}$
$f_{VL}(\phi _{V},\nabla \phi _{V},\phi _{L},\nabla \phi _{L})=\epsilon _{VL}(\phi _{L}\nabla \phi _{V}-\phi _{V}\nabla \phi _{L})^{2}+U_{VL}\phi _{V}^{2}\phi _{L}^{2}$

For simplicity, choosing solid-vpaor pair as an example, the variational derivative is,

${\frac {\delta F_{SV}}{\delta \phi _{S}}}=2U_{SV}\phi _{S}\phi _{V}(\phi _{S}-\phi _{V})-\nabla \cdot [\epsilon _{SV}(\mathbf {n_{S}} )|\phi _{V}\nabla \phi _{S}-\phi _{S}\nabla \phi _{V}|]$
       $-\sum _{i=1,2,3}{\frac {\partial }{\partial x_{i}}}\left[|\phi _{V}\nabla \phi _{S}-\phi _{S}\nabla \phi _{V}|^{2}\epsilon _{SV}(\mathbf {n_{S}} ){\frac {\partial \epsilon _{SV}(\mathbf {n_{S}} )}{\partial \phi _{,i}^{S}}}\right]$

where $x_{1}\equiv x$, $x_{2}\equiv y$, $x_{3}\equiv z$, and $\phi _{,i}^{S}\equiv \partial \phi _{S}/\partial x_{i}$.

The methods to obtain the explicit expression of $\epsilon _{SV}(\mathbf {n_{S}} )$ and ${\frac {\partial \epsilon _{SV}(\mathbf {n_{S}} )}{\partial \phi _{,i}^{S}}}$ were already discussed in the previous wiki page.

Algorithm

Similar to the previous wiki page,in matlab, We present liquid $\phi _{1}$, solid $\phi _{2}$, and vapor $\phi _{3}$, three phase fields by a $1\times 3$ cell. Each element phi{i}(m,n,p) of the cell is a three-dimensional array, of size $N_{x}\times N_{y}\times N_{z}$. i goes from 1 to 3, m goes from 1 to $N_{x}$, n goes from 1 to $N_{y}$, and p goes from $N_{z}$. For simplicity, a uniform grid size $h$ is used in all three directions., by default, periodic boundary conditions (PBC) are applied in all three directions.

For each phase $\phi _{i}$, the first derivatives, $\nabla _{x}\phi _{i}$, $\nabla _{y}\phi _{i}$, $\nabla _{z}\phi _{i}$, the second derivative, $\nabla ^{2}\phi _{i}$ and $(\nabla \phi _{i})^{2}$, are computed the following manner.

dphidx{i} = (phi{i}(pxind,:,:)-phi{i}(nxind,:,:))/(2*h);
dphidy{i} = (phi{i}(:,pyind,:)-phi{i}(:,nyind,:))/(2*h);
dphidz{i} = (phi{i}(:,:,pzind)-phi{i}(:,:,nzind))/(2*h);
d2phi{i}  = (phi{i}(pxind,:,:)+phi{i}(nxind,:,:) ...
             +phi{i}(:,pyind,:)+phi{i}(:,nyind,:) ...
             +phi{i}(:,:,pzind)+phi{i}(:,:,nzind)-6*phi{i})/(h^2);
dphi_SQR{i} = (dphidx{i}.*dphidx{i}) + (dphidy{i}.*dphidy{i}) + (dphidz{i}.*dphidz{i});

For the solid phase i=2,