View source for Single Phase Field Model with Anisotropic Surface Energy

<P ALIGN="CENTER">
<FONT SIZE="+3" color="darkred"><STRONG>
Single Phase Field Model with Anisotropic Surface Energy</STRONG></font></P>
<DIV>
<P ALIGN="CENTER"><STRONG>Seunghwa Ryu and Wei Cai</STRONG></P>
</DIV>

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). 
<HR>

==Download files==

The matlab program can be downloaded from this link: [[media:pf3d_single_aniso.m.txt‎ | 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/8/86/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 [[Single_Phase_Field_Model_with_Isotropic_Surface_Energy | the previous wiki page]], the fundamental degree of freedom in this model is a phase field, <math>\phi(\mathbf{x})</math>.  The only difference is that the gradient term in the the free energy density function now depends on the local orientation of the gradient.

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

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

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

 <math>\frac{\delta F}{\delta \phi} = 2 U \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] </math>

where <math>x_1\equiv x</math>, <math>x_2\equiv y</math>, <math>x_3\equiv z</math>, and <math>\phi_{,i}\equiv \partial \phi/\partial x_i</math>.

In this model, we use the following equation to describe the orientation dependence of <math>\epsilon</math> with cubic symmetry,

 <math>\epsilon(\mathbf{n}) = \epsilon_0 + \epsilon_1 \cdot (n_x^2n_y^2+n_y^2n_z^2+n_z^2n_x^2) + \epsilon_2 \cdot (n_x^2n_y^2n_z^2) + \epsilon_3 \cdot (n_x^2n_y^2+n_y^2n_z^2+n_z^2n_x^2)^2 </math>

where <math>\epsilon_0, \epsilon_1, \epsilon_2, \epsilon_3</math> are constants.  To obtain the derivative <math>\partial \epsilon(\mathbf{n})/\partial \phi_k</math>, we need to use the chain rule,

 <math>\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}</math>

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

 <math> n_i = \frac{\phi_{,i}}{\sqrt{(\phi_{,1})^2 + (\phi_{,2})^2 + (\phi_{,3})^2 }} </math>
 <math> \frac{\partial n_j}{\partial \phi_{,k}} = \frac{\delta_{jk}}{|\nabla \phi|} - \frac{\phi_{,j}\phi_{,k}}{|\nabla \phi|^3} </math>
 

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

==Algorithm==

Similar to [[Single_Phase_Field_Model_with_Isotropic_Surface_Energy | the previous wiki page]], we represent the field <math>\phi(\mathbf{x})</math> by a three-dimensional array, phi(i,j,k), of size <math>N_x \times N_y \times N_z</math>.  In Matlab, i goes from 1 to <math>N_x</math>, j goes from 1 to <math>N_y</math>, k goes from 1 to <math>N_z</math>.  For simplicity, a uniform grid size <math>h</math> is used in all three directions.  By default, periodic boundary conditions (PBC) are applied in all three directions.  The first derivatives, <math>\nabla_x \phi</math>, <math>\nabla_y \phi</math>, <math>\nabla_z \phi</math>, and the second derivative, <math>\nabla^2 \phi</math>, 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 <math>\mathbf{n}</math> 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 <math>\epsilon(\mathbf{n})</math> 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 <math>(\nabla \phi)^2 \epsilon(\mathbf{n}) \frac{\partial \epsilon(\mathbf{n})}{\partial \phi_{,i}}  </math>, 

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 <math>\phi(\mathbf{x}) = 0.5</math>, 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''.)