Single Phase Field Model with Anisotropic Surface Energy

Single Phase Field Model with Anisotropic Surface Energy

Seunghwa Ryu and Wei Cai

This tutorial describes the theory and numerical implementation of a 3D phase field model with anisotropic surface 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_single_aniso.m. For this 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/5/5b/Pf3d_single_aniso.m.txt -O pf3d_single_iso.m

and run it within Matlab.

To run the MD++ test cases, first you need to install and compile MD++ on your computer. Follow the instructions on [Single_Phase_Field_Model_with_Isotropic_Surface_Energy | the previous wiki page].

Theory

Similar to the previous wiki page, the fundamental degree of freedome in this model is a phase field, $\phi (\mathbf {x} )$. The only difference is that the gradient term in the the free energy density function now depends on the local orientation of the gradient.

$F[\phi (\mathbf {x} )]=\int f(\phi (\mathbf {x} ),\nabla \phi (\mathbf {x} ))d^{3}\mathbf {x}$
$f(\phi ,\nabla \phi )={\frac {1}{2}}(\epsilon (\mathbf {n} ))^{2}(\nabla \phi )^{2}+U\phi ^{2}(1-\phi )^{2}$

Notice that we have replaced the constant $\varepsilon$ in the previous formulation by ${\frac {1}{2}}(\epsilon (\mathbf {n} ))^{2}$, where $\mathbf {n}$ is a unit vector specifying the local orientation of $\nabla \phi (\mathbf {x} )$.

For this free energy density function, the variational derivative is,

${\frac {\delta F}{\delta \phi }}=2U\phi (1-\phi )(1-2\phi )-\nabla \cdot [\epsilon (\mathbf {n} )^{2}\nabla \phi ]-\sum _{i=1,2,3}{\frac {\partial }{\partial x_{i}}}\left[(\nabla \phi )^{2}\epsilon (\mathbf {n} ){\frac {\partial \epsilon (\mathbf {n} )}{\partial \phi _{,i}}}\right]$

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

In this model, we use the following equation to describe the orientation dependence of $\epsilon$ with cubic symmetry,

$\epsilon (\mathbf {n} )=\epsilon _{0}+\epsilon _{1}\cdot (n_{x}^{2}n_{y}^{2}+n_{y}^{2}n_{z}^{2}+n_{z}^{2}n_{x}^{2})+\epsilon _{2}\cdot (n_{x}^{2}n_{y}^{2}n_{z}^{2})+\epsilon _{3}\cdot (n_{x}^{2}n_{y}^{2}+n_{y}^{2}n_{z}^{2}+n_{z}^{2}n_{x}^{2})^{2}$

where $\epsilon _{0},\epsilon _{1},\epsilon _{2},\epsilon _{3}$ are constants. To obtain the derivative $\partial \epsilon (\mathbf {n} )/\partial \phi _{k}$, we need to use the chain rule,

${\frac {\partial \epsilon (\mathbf {n} )}{\partial \phi _{k}}}=\sum _{j}{\frac {\partial \epsilon (\mathbf {n} )}{\partial n_{j}}}\cdot {\frac {\partial n_{j}}{\partial \phi _{k}}}$

The following relations are needed to obtain the explicit expression of the variation derivative.

$n_{i}={\frac {\phi _{,i}}{\sqrt {(\phi _{,1})^{2}+(\phi _{,2})^{2}+(\phi _{,3})^{2}}}}$
${\frac {\partial n_{j}}{\partial \phi _{,k}}}={\frac {\delta _{jk}}{|\nabla \phi |}}-{\frac {\phi _{,j}\phi _{,k}}{|\nabla \phi |^{3}}}$

(the following are copied from pf3d_single_iso.m wiki page)

Algorithm

Similar to the previous wiki page, we represent the field $\phi (\mathbf {x} )$ by a three-dimensional array, phi(i,j,k), of size $N_{x}\times N_{y}\times N_{z}$. In Matlab, i goes from 1 to $N_{x}$, j goes from 1 to $N_{y}$, k goes from 1 to $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. The first derivatives, $\nabla _{x}\phi$, $\nabla _{y}\phi$, $\nabla _{z}\phi$, and the second derivative, $\nabla ^{2}\phi$, are computed by centered numerical difference.

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

We compute the local surface normal vector $\mathbf {n}$ as follows,

dphi_SQR = (dphidx.*dphidx) + (dphidy.*dphidy) + (dphidz.*dphidz);
absdphi = sqrt(dphi_SQR + 1e-24);
   
nx = dphidx ./ absdphi;
ny = dphidy ./ absdphi;
nz = dphidz ./ absdphi;

We then compute the local $\epsilon (\mathbf {n} )$ value, and the free energy functional,

triplesum = (nx.^2).*(ny.^2) + (ny.^2).*(nz.^2) + (nz.^2).*(nx.^2);
EPS_loc   = eps0 + eps1.*triplesum + eps2.*(nx.*ny.*nz).^2 ...
                    + eps3*(triplesum.^2);

F = (U * sum(sum(sum( phi.^2.*(1-phi).^2 ))) ...
    + 0.5 * sum(sum(sum( EPS_loc.^2 .* dphi_SQR)))) * (h^3);

The variational derivative is computed by first computing $(\nabla \phi )^{2}\epsilon (\mathbf {n} ){\frac {\partial \epsilon (\mathbf {n} )}{\partial \phi _{,i}}}$,

and then computing its spatial derivatives numerically,

   ...
   dEPS_dphidx = dEPS_dnx.*dnx_dphidx + dEPS_dny.*dny_dphidx + dEPS_dnz.*dnz_dphidx;
   dEPS_dphidy = dEPS_dnx.*dnx_dphidy + dEPS_dny.*dny_dphidy + dEPS_dnz.*dnz_dphidy;
   dEPS_dphidz = dEPS_dnx.*dnx_dphidz + dEPS_dny.*dny_dphidz + dEPS_dnz.*dnz_dphidz;
   ...
   tmp_x = (EPS_loc.^2).*dphidx + dphi_SQR.*EPS_loc.*dEPS_dphidx;
   tmp_y = (EPS_loc.^2).*dphidy + dphi_SQR.*EPS_loc.*dEPS_dphidy;
   tmp_z = (EPS_loc.^2).*dphidz + dphi_SQR.*EPS_loc.*dEPS_dphidz;    
   
   dFdphi = 2*U*phi.*(1-phi).*(1-2*phi) ...
          - (tmp_x(pxind,:,:)-tmp_x(nxind,:,:))/(2*h) ...
          - (tmp_y(:,pyind,:)-tmp_y(:,nyind,:))/(2*h) ...
          - (tmp_z(:,:,pzind)-tmp_z(:,:,nzind))/(2*h) ;

Again, the forward-Euler method is used for time integration for both the Ginzburg-Landau model and Cahn-Hillard model (depending on whether equation_type = 0, 1, or 2).

Test cases

We will focus on a test case in which an initially spherical interface become faceted as it evolves with time. This can be accomplished by typing the following commands in Matlab,

clear all
equation_type = 2; dt = 1e-3; pf3d_single_aniso

Here equation_type = 2 specifies the use of Ginzburg-Landau equation but with the total amount of solid phase constrained (through a Lagrangian multiplier). dt = 1e-3 specifies the time step. You should see the isosurfaces where $\phi (\mathbf {x} )=0.5$, depicting the solid-liquid interface, which evolves with time. Typing pf3d_single_iso again (without clear all) continues the simulation for another Niter = 5000 steps.

(MD++ test cases. The code is ready. Just need to prepare the tcl script and write up the description. In MD++, we cannot plot isosurfaces, but can only visualize the phase field by plotting a circle at grid points where the phase field value lies within a window (e.g. [0.45, 0.55]). However, MD++ can save simulation snapshots into .cn files, which can be loaded into matlab and plotted using plot_mdpp_phasefield.m.)